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Abstract 
A  FORTRAN  program  is  presented  that  corrects  for  collimation 
effects  in  small  angle  x-ray  scattering  measurements.   The  method 
employed  in  calculations  is  based  on  numerical  solution  of  an 
integral  equation,  which  is  written  as  a  Volterra  equation  of  the 
first  kind.   This  equation  is  reduced  to  a  system  of  simultaneous 
equations,  which  can  be  solved  by  back  substitution.   Several 
different  algorithms  are  tried,  based  on  different  interpolating 
polynomials  for  I (x) ,  the  true  scattering  intensity,  which  is 
expressed  as  a  function  of  the  scattering  angle  x.   The  problem  of 
numerical  stability,  which  is  inherent  to  the  computational  algo- 
rithms employed  in  the  algebraic  solution  equation  is  discussed  and 
exemplified  with  actual  computations. 

Keywords:   Desmearing;  FORTRAN  program;  integral  equation;  inter- 
polating polynomial;  numerical  stability;  slit  correction;  small 
angle  scattering. 


CHAPTER  1 
1 .   INTRODUCTION 

Measurement  of  small  angle  scattering  of  x-rays  is  an  extremely 
valuable  tool  for  the  study  of  macromolecular  systems.   In  recent  years, 
these  measurements  have  been  extensively  employed  in  order  to  determine 
the  three  dimensional  structure  of  triblock  copolymers.   These  coploymers 
form  domains  whose  size  is  ideally  suited  for  small  angle  x-ray  scatter- 
ing studies.   As  a  consequence,  the  composition  of  the  domains  and  the 
changes  of  their  morphology  with  variation  of  stress,  temperature  or 
pressure  can  be  easily  studied.   Since  triblock  copolymers  are  fre- 
quently used  in  the  plastics,  rubber  and  coating  industry,  their  study 
and  characterizations  by  small  angle  scattering  is  of  prime  technologi- 
cal importance. 

In  x-ray  scattering  measurements,  narrow  slits  are  often  used  in 
order  to  collimate  the  incident  beam  sufficiently  to  allow  scattering 
measurements  close  to  the  incident  beam.   Compared  to  a  pinhole  of 
diameter  approximately  equal  to  the  slit  width,  the  slit  system  increases 
the  available  scattering  energy.   Therefore,  the  experimental  measure- 
ments of  the  angular  distribution  of  the  intensity  of  small  angle  x-ray 
scattering,  instead  of  giving  the  scattering  at  a  single  angle,  represent 
the  average  intensity  over  a  range  of  angles.   Because  of  this,  the 
measured  scattering  curve  is  often  very  different  from  the  curve  ob- 
tained with  point  collimation.   Hence,  the  collimation  effect  results  in 
"smearing"  of  the  scattered  radiation. 


When  the  slit  is  very  narrow,  the  effect  of  the  width  of  the  slit 
on  the  smearing  of  x-ray  scattering  can  be  neglected,  and  only  the 
effect  of  the  slit  length  needs  to  be  considered.   In  some  cases,  the 
slit  width  produces  changes  in  the  scattering  curves  that  should  not  be 
neglected.   Since  the  distortion  caused  by  the  width  of  the  collimating 
slit  is  small,  the  correction  for  the  width  effects  might  be  well  approxi- 
mated using  the  method  developed  by  Taylor  and  Schmidt  [1].   In  this 
method  the  measured  scattering  curve  is  corrected  for  the  width  effect, 
before  it  is  corrected  for  the  length  effect.   In  the  present  work,  the 
method  of  reference  [1]  is  employed.   However,  more  than  one  width  slit 
correction  is  being  applied  to  the  measured  scattering  curve.   This  is 
done  in  order  to  extend  the  method  of  reference  [1]  to  cases  when  the 
width  of  the  collimating  slits  is  altered  at  more  than  one  scattering 
angle. 

The  extent  of  the  angular  region  over  which  the  collimation  system 
averages  or  "smears"  the  scattered  radiation  can  be  expressed  by  a 
weighting  function.   This  weighting  function  represents  the  intensity 
distribution  within  the  primary  beam  in  the  longitudinal  direction.   In 
many  cases,  with  the  counter  tube  measurements,  the  length  of  the  slit 
with  which  the  intensity  distribution  is  measured  also  has  an  influence 
on  the  scattered  intensity.   If  the  slit  length  is  equal  to  k,  then  the 
corrected  weighting  function  is  equal  to  the  convolution  of  the  un- 
corrected weighting  function  with  a  step  function  which  represents  the 
length  of  the  slit.   Thus,  the  corrected  weighting  function  W  is  related 
to  the  uncorrected  weighting  function  w  through 
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(t)  =   /w(u)h(t-u)du 

—00 

(1) 


h(u)  =1       0$  u  <k 

h(u)  =0       k<  u  <°° 
w  is  a  complicated  function  of  the  geometrical  parameters  of  the 
collimating  system.   An  excellent  review  on  the  method  to  calculate  the 
weighting  function  from  the  geometry  of  the  system  is  presented  by 
Hendricks  and  Schmidt  [2],   There  is  another  way  to  determine  the 
weighting  function  with  reasonable  accuracy  without  having  to  compute  it 
from  the  geometry  of  the  collimating  system.   In  this  method,  the 
weighting  function  is  determined  by  measuring  the  intensity  of  the  beam 
in  the  plane  of  registration.   This  is  done  by  taking  various  intervals 
equal  to  the  length  of  the  scanning  slit  at  the  detector. 

The  integral  equation  which  relates  the  experimentally  determined 
intensity  I (x) ,  measured  at  the  scattering  angle  x,  to  the  true  scatter- 
ing intensity  I (x)  and  to  the  weighting  function  W(y)  is  [3]: 

oo 

I(x)  -  /~W(y)  I((x2  +  y2)1/2)dy.  (2) 
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There  is  a  comprehensive  literature  which  deals  with  solving  integral 
equation  (2)  [4] .   Most  of  the  methods  reported  involve  making  certain 
assumptions  about  the  functional  form  of  W(y)  and  I(x).   An  exception  is 
the  procedure  due  to  Vonk  [5]. 

In  reference  [6]  the  author  proposed  a  class  of  multistep  methods 
which  like  Vonk's  does  not  require  assumptions  about  the  functional  form 


of  I  and  W.   These  methods  start  with  eq.  (2)  rewritten  for  positive 
scattering  angle  x  as  a  linear  Volterra  integral  equation  of  the  first 
kind,  [6] 


h(u)  =  /*K(u,x)I(x)dx  0<u$x<°° 
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The  experimentalist  provides  smeared  intensity  data  I(xq),  I(xi), 

.  .  .  I(x  )  where  0<a  =  x„  <  xi  .  .  .  <x  =  b  <  °°.   Generally,  the  x. 
n  0     1         n  J  '      l 

are  equally  spaced,   A  square  grid  is  set  up  on  [a,b]  x  [a,b]  of  the 
u-x  plane  with  mesh  size  Ax  and  mesh  points  {(u.,  x.)  =  (iAx,  jAx), 
m£i,j<:N}.   mAx  is  the  lowest  scattering  angle  at  which  an  intensity 

'Xj 

measurement  is  taken.   Interpolation  is  used  to  obtain  values  of  I(xn) 
if  measured  data  at  angle  xn  are  not  available.   Quadrature  rules  are 
then  used  to  discretize  eq.  (3).   Accordingly, 
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h(iAx)    /  j  a. .  K(iAx,  JAx)  I  (jAx) 
j=i  1J 


(4) 


i,j  =  m,m  +  1,  .  .  .  N-l. 


To  complete  the  set  of  linear  equations  we  proposed  in  reference  [6] 
to  use  the  equation 


h'(NAx)  =  j   [h(NAx)-h((N-l)Ax)]=K(NAx,  NAx)I(NAx) .    (4a) 

Li  A 

Here  I(jAx)  denotes  an  approximation  to  I(jAx).   The  set  of  values 
I(jAx),  j=m,  .  .  .  ,  N  will  be  treated  as  components  of  the  solution 
vector  I.   Ways  to  approximate  I(jAx)  and  choices  of  quadrature  rules 
are  discussed  in  Chapter  3  while  the  evaluation  of  h  is  discussed  in 
Chapter  4,  Section  E.   For  simplicity,  we  will  frequently  drop  the 
bar  from  I(jAx)  when  denoting  approximate  values. 

The  step-by-step  discretization  of  eq.  (3)  leads  to  a  convergent 
method  if  in  rather  loose  terms  the  difference  between  the  approximate 
value  I(iAx)  and  the  true  value  I(iAx)  tends  to  zero  as  Ax-K)  with 
x=iAx  fixed  and  xe[a,b]. 

Integral  equations  of  the  first  kind  are  known  to  be  ill-posed  [7]. 
This  means  that  small  perturbations  in  h(u)  or  in  K(u,x)  may  have  a 
large  effect  on  the  solution  \.      Translated  into  our  physical  problem, 
one  could  suspect  that  small  errors  in  the  experimental  measurements  of 
I(x)  or  in  the  collimation  weighting  function  W  will  be  greatly  enhanced 
in  the  unsmearing  process.   The  extreme  sensitivity  of  the  solution 
vector  1  to  errors  in  the  experimental  data  can  be  greatly  reduced  by 
converting  the  Volterra  integral  equation  of  first  kind  into  the  equa- 
tion of  the  second  kind  by  differentiation.   The  disadvantage  in  making 
such  a  conversion  is  that  numerical  differentiation  of  experimental  data 
can  introduce  serious  errors.   The  method  for  desmearing  experimental 
data  presented  in  the  Technical  Note  should  be  employed  in  cases  when 
numerical  differentiation  of  experimental  data  is  to  be  avoided. 


There  are  several  sources  for  possible  numerical  errors.   The  first 
one  is  inherent  in  the  quadrature  rule  employed  to  discretize  the  integ- 
ral equation.   The  second  one  is  random  errors  or  noise  present  in 
experimental  data.   Actually,  these  two  sources  of  errors  cannot  be 
treated  separately.   Certain  quadrature  rules  (e.g.  trapezoid  rule) 
require  that  the  data  be  "sufficiently  smooth"  to  insure  numerical 
stability.   (A  function  is  defined  to  be  "sufficiently  smooth"  if  it  is 
continuously  dif f erentiable  up  to  sufficiently  high  order  with  respect 
to  its  argument  [8].)   However,  certain  quadrature  methods  will  always 
be  unstable.   To  get  an  insight  into  numerical  instability,  assume  that 
eqs.  (4-4a)  are  solved  by  back  substitution,  starting  with  I(NAx).   The 
method  is  absolutely  unstable  if  the  successive  values  of  I(x)  tend  to 
diverge  as  Ax  is  made  smaller,  while  keeping  NAx  constant.   In  this  work 
several  algorithms  for  discretization  of  the  integral  equation  will  be 
employed  and  which  are  not  absolutely  unstable.   (These  algorithms  along 
with  a  stability  criterion  will  be  presented  in  Chapter  3.)   The  purpose 
of  presenting  several  algorithms  is  to  enable  the  user  of  the  program  to 
compare  results  based  on  various  algorithms  and  to  determine  whether  the 
maxima  and  minima  present  in  the  unsmeared  scattering  curve  are  due  to 
random  errors  or  whether  they  are  real. 

The  program  presented  in  this  Note  does  not  provide  for  the  smooth- 
ing of  the  experimental  data.   This  should  be  the  user's  responsibility. 
In  particular,  the  use  of  the  trapezoidal  rule  has  been  known  to  lead  to 
even-odd  oscillations  in  the  computed  discrete  set  of  values  for  I(x) 
[8],   These  oscillations  can,  however,  be  reduced  by  smoothing  the 
results.   The  method  for  smoothing  the  computed  I(x)  employed  in  this 
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work  is  a  modification  of  the  method  suggested  by  Linz  [9] ,  and  is 
described  in  Chapter  3.   We  also  employ  an  algorithm  for  discretization 
of  the  integral  equation  which  is  based  on  forward  differences  in  a  two- 
step  quadrature  method  and  designated  as  the  Nystrom  algorithm.   We 
found  that  with  this  algorithm,  no  smoothing  of  the  results  was  neces- 
sary, since  the  computed  values  for  I(x)  were  relatively  free  of  even- 
odd  oscillations. 

There  is  another  source  of  error  that  is  brought  about  by  the 
termination  of  the  calculations  at  a  certain  point.   This  termination 
introduces  a  finite  upper  limit  to  the  integral  equation  (3) .   The 
perfect  collimation  theory  calls  for  an  infinite  upper  limit  in  this 
integral  equation.   This  error  will  be  more  noticeable  in  the  outer  part 
of  the  scattering  curve,  close  to  its  terminal.   Therefore,  for  the 
results  to  be  physically  meaningful,  it  is  important  to  extend  the 
measurements  as  far  as  possible  beyond  the  significant  range  of  the 
scattering  angles.   We  did  find  that,  if  this  precaution  is  taken,  the 
termination  error  will  not  affect  the  locations  and  the  relative  inten- 
sities of  the  maxima  and  minima  in  the  desmeared  intensity  curve. 

To  summarize,  the  procedure  employed  in  this  work  for  the  data 
desmearing,  based  on  a  class  of  step-by-step  finite  difference  methods, 
has  certain  advantages  and  disadvantages  over  the  procedures  reported  in 
the  pertinent  literature.   The  advantages  are  two-fold:   1)   No  numeri- 
cal differentiation  of  the  experimental  data,  a  process  which  may  intro- 
duce serious  error,  is  required;  2)   No  restriction  on  the  functional 
form  for  the  weighting  function  is  imposed.   The  weighting  function  can 

be  presented  as  a  set  of  data  which  are  predetermined  by  the  experiment 
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either  by  intensity  measurements  in  the  plane  of  registration  or  by 
calculation  from  the  geometry  of  the  collimating  system. 

The  main  disadvantage  of  the  method  of  solving  eq.  (3)  by  finite 
differences  is  that,  even  with  sufficiently  smooth  data,  many  quadrature 
methods,  commonly  employed  in  solving  finite  difference  equations,  will 
lead  to  divergent  results  as  Ax-K). 

In  this  work  we  provide  a  few  tests  of  the  method  of  solving  eq. 
(3).   These  tests  are  of  two  kinds:   1)   Solution  checking  and;  2) 
testing  the  various  algorithms  against  test  functions  for  I(x),  which 
lead  to  analytically  known  perfect  collimation  intensities  I(x).   The 
solution  checking  calculations  are  incorporated  in  the  program.   Two 
test  functions  for  I(x)  are  employed.   One  is  a  Gaussian  function.   The 
other  represents  a  scattering  curve  for  spheres  with  uniform  electron 
densities.   This  function  has  been  tabulated  by  Schmidt  [10].   For  the 
first  test  function,  the  perfect  collimation  intensities  are  also  of 
Gaussian  form.   The  second  test  function  corresponds  to  the  perfect 
collimation  intensity  of 

sin(xa)  -  xacos(xa) 


I(x)  =  A 


(xa) 


where  a  is  the  radius  of  the  scattering  sphere  and  x  is  proportional  to 
the  scattering  angle.   In  both  test  functions,  a  constant  weight  func- 
tion was  assumed.   Both  test  functions  lead  to  satisfactory  results. 

The  theoretical  background  for  the  methods  employed  in  conjunction 
with  the  numerical  solution  of  the  set  of  simultaneous  equations  (eqs. 
4-4a)  and  the  ways  to  cope  with  the  problem  of  numerical  stabilities  is 


reviewed  in  Chapter  3.   In  Chapter  4,  the  actual  computations  are  des- 
cribed, followed  by  the  FORTRAN  program  (Chapter  5).   In  Chapter  6,  the 
various  algorithms  employed  in  construction  of  numerically  convergent 
methods  and  their  effect  on  the  solution  vector  I  are  discussed  and 
shown  in  Figures  1-4.   The  experimental  data  for  this  investigation  were 
provided  by  J.  J.  Weeks  of  this  laboratory.   In  Figure  5  we  show  how 
unsmearing  of  data  lead  to  a  formation  of  peaks  in  the  true  intensity 
curve.   The  data  for  this  Figure  were  supplied  by  D.  Mclntyre  of  Polymer 
Institute,  University  of  Akron.   At  the  end  of  Chapter  6,  calculations 
based  on  the  two  test  functions  described  above,  are  shown  (Table  5  and 
Figure  6) . 

For  the  person  who  is  interested  only  in  the  preparation  and  use  of 
the  program,  we  present  in  Chapter  2  the  organization  of  the  FORTRAN 
program,  and  provide  the  instructions  for  preparation  of  the  input  data 
required  for  the  actual  computations. 
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CHAPTER  2 
1.   PROGRAM  ORGANIZATION 
A.   Available  Programs 

There  are  two  separate  programs.   The  first  program,  listed  in 
Chapter  5,  should  be  used  when  N,  the  total  number  of  experimental  data 
employed  in  calculations  exceeds  150.   In  this  case,  storing  matrices 
would  require  very  large  computer  storage.   To  avoid  the  storage  problem, 
the  matrix  elements  are  computed  as  functions  for  each  pair  of  indices  i 
and  j,  with  i<j<:N.   This  process  is  more  time-consuming,  than  with  the 
matrix  storage.   Therefore,  for  N>150,  if  computer  time  is  of  prime 
importance,  the  weighting  function  should  be  approximated  by  a  trapezoid 
or  by  a  constant,  if  possible.   The  case  where  the  weighting  function  is 
represented  by  a  set  of  data  is,  by  far,  the  most  time-consuming. 

The  second  program  is  designed  for  N<150.   Since  the  components  of 
the  matrix  employed  In  calculations  (see  Chapter  4  for  details  on  the 
matrices  employed)  are  stored  in  the  computer  memory  at  the  beginning  of 
the  computations,  they  can  be  used  in  conjunction  with  various  algorithms 
employed  in  the  computation  of  the  solution  vector  1.   Except  for  this 
detail,  the  two  programs  are  identical.   This  second  program  is  available 
upon  request.   It  is  not  listed  here. 

B.   Instructions  and  Input  Data 

The  program  requires  that  the  components  of  the  input  vector  1 


a- 

[ 

o 

be  spaced  at  equal  intervals  of  x.   This  increment,  Ax,  is  customarily 


expressed  in  milliradians.   Very  often,  the  experimentalist,  at  some 
point,  increases  the  spacings  between  two  consecutive  scattering  angles, 
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Suppose  that  at  a  certain  scattering  angle  equal  to  pAx,  the  spacing  be- 
tween scattering  angles  is  increased  from  Ax  to  £Ax,  where  H  is  a  integer. 
The  computer  program  will  interpolate  by  quadratic  forward  interpolation 
the  input  data,  starting  with  the  angle  pAx  to  give  I  values  at  angular 
increment  Ax.   The  last  interval,  between  x=(N-l)Ax  and  x=NAx,  is  inter- 
polated linearly.   The  new  input  vector  1  will  contain  N  =N  +  (N  -p) 
f  J  r  too 

(£-1)  equally  spaced  angles,  where  N   is  the  number  of  input  values  fed 
into  the  machine.   Only  one  change  in  data  spacing  is  allowed. 

For  example,  suppose  that  the  measured  data  for  the  scattering 
intensity  were  initially  collected  at  a  constant  interval  of  Ax=0.2 
mrad.   Suppose  that  the  experimentalist,  starting  with  a  scattering 
angle  equal  to  20  mrads ,  increased  the  spacing  between  consecutive 
measurements  from  0.2  mrads  to  1  mrad.   If  the  initial  number  of  ex- 
perimental measurements  was  200,  N  will  be  equal  to  200+(200-100) • (4)=600 , 
since  p=100  and  1=5.      The  largest  scattering  angle  is  equal  to  120 
mrads. 

The  lowest  scattering  angle,  at  which  the  initial  measurement  is 
taken,  ordinarily  exceeds  the  spacing  Ax  between  two  consecutive  scatter- 
ing angles.   Let  m  be  the  ratio  of  the  lowest  scattering  angle  to  the 
spacing  between  two  consecutive  scattering  angles,  rounded  to  the  nearest 
integer.   As  a  consequence  of  this  rounding  off,  the  measurements  of 
scattering  angles  might  be  off  by  a  much  as  Ax/2.   To  avoid  this  error, 
the  lowest  scattering  angle  should  be  chosen  in  such  a  way  that  it 
equals  an  integer  number  of  the  spacing  between  scattering  angles. 
Thus,  mAx  measures  the  lowest  scattering  angle  from  which  the  measure- 
ments  are  taken.   The  array  which  represents  the  input  vector  ]  after 

12 


it  is  read  into  the  computer,  is  shifted  by  m  locations. 

Let  there  be  gi  values  at  initial  interval  Axj  and  g2  values  at 
augmented  interval  Ax2 ,  and  let  £=Ax2/Ax]_.   The  program  will  fill  vector 

and  generate  vector  X  so  that  one  interval  Ax  will  apply  everywhere. 
Only  one  interpolation  is  allowed.   We  have: 

g!  =  p-m+1,  g2=No-gi . 
The  following  diagram  1  (Diagram  1  is  presented  on  page  I)  shows  the 
preparation  of  the  augmented  input  vector  T .   In  diagram  1,  x  =mAx, 
x =(m+l)Ax,  etc.   Starting  angle  is  equal  to  mAx. 


'Xj 


Card  1  provides  four  input  data  needed  to  fill  vector  T  .   These 
are:  1)   N  ,  total  number  of  input  data  (=g]+g2) >  2)   m,  defined  so  that 
mAx  measures  the  lowest  scattering  angle  from  which  measurements,  taken 
at  equal  interval  of  Ax,  are  taken;  3)   p,  defined  so  that  pAx  is  the 
starting  angle  for  change  in  data  spacing  and;  4)   the  number  of  inter- 
polated values  required  in  each  Ax2  interval.   This  number  is  equal 
to  1-1,   where  £=Ax2/Ax1. 

These  four  input  parameters  are  read  as  NN,  ISTART,  IENTR  and  NSUB. 
They  are  punched  on  a  single  card,  FORMAT  4110. 

If  ISTART=0,  1}  is  the  intensity  at  the  zero  scattering  angle 
(might  be  obtained  by  extrapolation  of  data  to  zero  angle) . 

The  experimentalist  might  increase  the  width  of  the  slit  at  certain 

scattering  angles,  in  order  to  provide  a  sufficient  scattering  intensity 

at  larger  angles.   In  the  computer  program,  if  the  card  which  follows 

the  first  card  is  blank,  no  correction  for  the  slit-width  is  performed. 

For  each  slit-width  correction,  a  single  card  is  read  which  provides  two 

sets  of  data:   the  lower  limit  of  the  scattering  angle  (in  Ax  units) , 
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from  which  the  slit-width  correction  applies,  and  the  value  of  M2 ,  the 
second  moment  of  the  slit-width  weighting  function  W  ,  defined  as: 


w 


M2  =  f  t2W  (t)dt 
J  W 


W  is  normalized,  so  that 
w 


c 


W  (t)dt=l, 
w 


M2  (customarily  given  in  mrads2)  has  to  be  given  in  (Ax)z  units. 

A  separate  subroutine  performs  a  slit-width  correction  on  the  data 
after  the  interpolation.   The  correction  is  made  according  to  the  equation 

I(x)=I(x)-  |M2I"(x) 

Numerical  evaluation  of  the  second  derivative  of  I(x),  based  on  a  least- 
squared  smoothing  technique,  is  given  in  reference  [1].   Let  us  assume 
that  a  series  of  k  slit-width  corrections  apply,  each  in  a  different 
range  of  scattering  angles.   Let  the  second  moment  of  the  slit-width 
weighting  functions  in  the  range  WjAx^x^W^Ax  be  M2(l),  in  the  range  W2Ax< 
x<W3Ax  be  M2(2),  etc.   The  last  correction  always  applies  in  the  range 
W,  Ax£x<NAx.   The  following  diagram  2  (Diagram  2   is  presented  on  Page  I) 
illustrates  these  corrections. 

There  are  k  cards  for  slit-width  correction,  one  card  for  each 
range  of  scattering  angles  for  which  M2(i),  i=l,  2,  ...  k,  is  constant. 
These  cards  follow  the  first  card  of  the  program.   These  k  cards  are 
followed  by  a  blank  card.   A  blank  card  must  be  inserted,  whether  the 
data  are  corrected  for  the  slit-width  effect  or  not.   Each  one  of  these 
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k  cards  contain  two  numbers:   W.  and  M2(j),  j=l,  2,  .  .  . ,  k.   They  are 
punched  FORMAT  110,  F10.5. 

A  word  of  warning:   It  is  advisable  not  to  correct  data  for  the 
slit-width  effect  in  the  steep  part  of  the  scattering  curve,  since  the 
second  derivative  might  be  very  large.   It  is  recommended  to  apply  the 
slit-width  correction  only  at  scattering  angles  in  excess  of  5  milli- 
radians. 

The  program  user  has  four  directive  parameters  which  control  the 
computations  for  the  slit- length  effects.   These  four  parameters  are 
punched  on  a  single  card,  FORMAT  4110.   These  are:   MD,  IDC,  INV,  MM. 
MD  is  an  integer  which  must  have  one  of  the  values  1,  2,  3,  or  4.   This 
integer  controls  the  form  of  the  weighting  function  to  be  employed  in 
computations.   These  are,  in  order  indicated  by  the  value  given  to  MD: 
1)   Constant  weighting  function;  2)   Gaussian  weighting  function;  3) 
Data  weighting  function  and;  4)   Trapezoid  weighting  function.   Accord- 
ing to  the  form  of  the  weighting  function,  additional  information  must 
be  supplied,  which  describes  the  shape  and  form  of  this  function,  as 
follows : 

MD=1  No  additional  input  data  are  needed 
MD=2   The  weight  function  is  defined  as 

W.  =  exp  -  [(iA)2/G] 
A  is  the  interval  between  successive  values  of  W  and  G/2  is  its  variance, 
A  is  determined  in  such  a  way  so  that  a  table  of  21  values  of  W  are 
computed, 

W.=exp-(0.01i2);  i=0,  1,  ...,  20. 
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Thus,  A,  the  interval  between  two  successive  values  of  W,  is  equal  to  0.1 

/G. 

-4 
W  is  cut-off  when  it  is  less  than  e   .   A  single  card  giving  G, 

FORMAT  F10.2,  must  be  inserted. 

MD=3 :   The  data  weighting  function  is  given  by  a  set  of  values  W(R) .   R 

should  be  equally  spaced,  but  the  spacing  can  be  larger  than  the  spacing 

of  x  in  I(x).   R  is  given  in  Ax  units.   W  can  be  given  in  arbitrary 

units.   A  set  of  cards  giving  the  pairs  of  R  and  W,  one  pair  per  card, 

must  be  inserted,  FORMAT  F9.3,  F11.5.   As  many  cards  as  needed  can  be 

inserted.   Last  card  must  be  blank. 

MD=4:   A  single  card  giving  S}  and  S2  must  be  inserted  (FORMAT  2F10.2). 

Si  must  be  smaller  than  S2.   (The  trapezoid  is  defined  so  that  W=l  for 

r<Si  and  W=0  for  r<S2,  and  it  is  symmetric  around  r=0)  (illustrated  on 

page  33).  S\    and  S2  are  in  Ax  units. 

IDC,  which,  if  different  from  zero,  causes  the  computer  to  print  out  the 

values  of  the  integrals 

NAx 


1  '     J      (x2-(iAx)2)V2 


iAx 

for  all  values  of  i  from  m  to  N. 

INV,  which,  if  different  from  zero,  causes  the  input  vector  f  to  be 
recomputed  from  the  output  vector    (inversion  of  calculations) . 
MM,  which,  if  different  from  zero,  causes  the  integrals  h.  to  be  cal- 
culated by  the  lower-order  approximation  (see  Chapter  4  for  details). 
The  input  vector  t  is  read  after  all  the  above  instructions  were 

supplied.   For  MD=1  the  input  vector  follows  immediately  the  blank 
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card,  which  terminates  the  information  needed  for  the  slit-width  correc- 
tion. 

For  MD^l,  the  input  vector  follows  cards  which  characterize  the 
particular  weighting  function  employed  in  calculations. 

The  input  vector  ''  is  represented  by  a  series  of  real  numbers, 
eight  numbers  per  card,  FORMAT  8F10.0.   The  components  of  the  input 

'Xj  Oi     'Xj  'Xj  'Xj 

vector  T  are  Ij,  I2,  •  ..,  IN  •   After  N  values  of  T  are  fed  into  the 

o 
machine,  the  computations  will  start,  even  if  there  are  more  data 

cards  to  be  read.   At  the  conclusion  of  the  computations,  new  sets  of 

data  cards  will  be  read,  and  new  computations  will  be  initiated.   The 

program  can  handle  more  than  one  set  of  data,  however,  the  instructions 

and  specifications  must  remain  the  same.   The  computations  are  stopped 

'Xj 

when  a  blank  card  is  read  for  Ij. 

The  instructions  for  the  input  parameters  and  variables  are 
summarized  in  Table  1. 

C.   Brief  Description  of  the  Main  Program  and  its  Subroutines 
The  main  program  reads  all  the  necessary  parameters  and  variables 
which  control  the  computation,  followed  by  the  input  data.   The  input 
data  are  first  shifted  and  interpolated  as  needed.   The  main  program 
calls  subroutine  SWIDTH  for  the  corrections  required  for  the  effects  of 
the  variable  widths  of  the  collimating  slits.   The  calculations  performed 
by  the  subroutine  follow  the  method  of  Taylor  and  Schmidt  [1],   Next, 
subroutines  CGEN  and  BGEN  are  called  to  calculate  the  coefficients 
y.  ,  and  b.    in  the  set  of  linear  equations 

N      ^      N 

£    y..i.  =    £    b..i. 
j=i     1J  J       j=i     XJ  J 
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or,  in  matrix  form,  rT=DT .   The  main  program  performs  an  inversion  of 
the  matrix  equation  to  compute  T .   The  main  program  also  prints  out  the 
results  and  performs  the  back-check  (optional) .   Subroutine  BGEN  calls 
subroutine  WGHT23  if  MD=2  or  3,  or  calls  subroutine  WHGT4  if  MD=4.   The 
subroutines  WGHT23  and  WGHT4  calculate  K(u,x) ,  the  kernel  of  the  integral 
eq.  (3) ,  employing  the  selected  weighting  function.   These  calculations 
are  performed  for  each  pair  of  indices  (i,j)  (in  notation  of  eq.  (3), 
u=iAx,  x=jAx,  and  m<(i,j)<N,  where  mAx  is  the  lowest  angle  at  which 
measurements  were  taken) . 

In  addition,  the  program  includes  two  function  subroutines,  DGEN 
and  CGEN.   These  subroutines  transform  the  coefficients  b.  .,  calculated 
for  the  rectangle  rule  approximation  (for  which  I(x)  is  constant  in  each 
integral  and  interpolates  one  end  point)  into  the  coefficients  which  ob- 
tain when  I(x)  is  replaced  by  certain  interpolating  polynomials.   This 
procedure  is  explained  in  detail  in  Chapters  3  and  4. 

D.   Comments  on  Units  Employed 

%      .  . 
T ,  T ,  and  ^  are  given  in  arbitrary  units  (usually  as  counts  per 

a, 
seconds).   The  coordinates  of  ]  and  of  U  are  given  in  units  of  Ax,  where 


Ax=(0. , n-Q.)  and  9.  is  the  scattering  angle  at  the  point  i.  If  Ax  is 
equal  to  s  mrads,  all  other  linear  coordinates  must  be  given  in  terms  of 
multiples  of  s.  No  explicit  value  for  Ax  is  needed,  since  in  machine 
calculations,  Ax=l.  The  variance  in  the  gaussian  distribution  and  M2 , 
the  second  moment  of  the  width  weighting  function  are  in  (Ax)2  units. 
For  example,  if  M2  is  measured  as  0.15  mrads2,  and  Ax  is  equal  to  0.2 
mrads,  then  M2 ,  in  (Ax)2  units,  is  equal  to  3.75. 
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Occasionally,  the  unsmearing  of  data  might  be  restricted  to  larger 
scattering  angles,  or  the  calculations  for  larger  scattering  angles  are 
to  be  performed  separately  from  the  calculations  at  very  low  scattering 
angles.   Suppose  that  the  calculations  start  with  the  scattering  angles 
of  pAx,  at  which  the  spacing  between  the  consecutive  scattering  angles 
is  equal  to  £Ax.   In  this  case,  ISTART=p,  and  no  interpolation  of  input 
data  is  needed.   All  the  linear  coordinates  employed  will  be  in  units  of 
£Ax  instead  of  in  units  of  Ax.   In  machine  calculations,  £Ax=l.   Simi- 
larly, M2,  the  second  moment  of  the  slit-width  weighting  function,  will 
be  in  (£Ax)2  units. 

If  the  scattering  intensities,  starting  at  certain  angle  pAx,  are 
much  smaller  than  the  intensities  at  the  lowest  measured  angles,  then 
the  unsmearing  of  data  for  very  small  scattering  angles  (0<pAx)  might  be. 
performed  on  a  restricted  range  of  the  scattering  angles,  terminated  at 
the  angle  pAx.   This  artificial  cutting  off  the  calculations  generally 
will  not  affect  the  pattern  of  the  scattering  curve,  that  is  the  rela- 
tive positions  and  intensities  at  various  maxima  and  minima.   In  this 
case,  the  computations  should  be  performed  separately  for  very  low 
scattering  angles  (9<pAx)  and  separately  for  larger  scattering  angles 
(0>pAx) .   This  procedure  might  be  preferable  if  the  computer  time  is  of 
primary  concern.   By  performing  calculations  separately  for  different 
ranges  of  the  scattering  angles,  N,  the  total  number  of  the  input  data, 
is  reduced.   Since  the  processor  time  involved  is  roughly  proportional 
to  N2 ,  considerable  computer  time  can  be  saved. 
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CHAPTER  3 
1.   THEORETICAL  BACKGROUND 

The  equation  which  relates  the  experimentally  determined  intensity 
of  scattered  rays    to  the  true  scattering  intensity  ','  is  given  by  eq. 
(2).   In  eq.  (2),  W  may  be  given  either  in  analytic  form  or  as  a  set  of 
discrete  values.   Eq.  (2)  is  transformed  as  follows:   both  sides  of  eq. 
(2)  are  multiplied  by  x/(x2-u2)1'2  and  integrated  over  x  from  u  to  °°. 
After  exchanging  the  orders  of  integration  one  obtains  eq.  (3).   The 
total  number  of  experimental  data  for  I(x)  is  N,  which  are  given  at 
equally  spaced  intervals  Ax.   We  assume  that,  if  I(x)  is  sufficiently 
small  when  x>NAx,  so  is  I(x).   Therefore,  both  the  upper  limits  of 
integrals  in  eq.  (3)  are  replaced  by  NAx.   Both  I(x)  and  I(x)  represent 
a  discrete  set  of  function  values  I.  and  I.,  respectively,  with  i  =  0, 
1,  ...,  N.   Thus,  i=0  corresponds  to  u=0,  i=N  corresponds  to  u=NAx. 
Actually,  the  calculations  seldom  start  with  u=0  but  from  a  small  angle 
equal  to  mAx,  m>0,  since  l(x=0)  is  not  known,  unless  it  is  obtained  by 
extrapolation  of  data  to  zero  angle. 

The  l.h.s.  of  eq .  (3)  can  be  written  as  a  :>et  of  N  linear  equa- 
tions: 

NAx  N 


h±  =   /   I(x)G(u,x)dx^  V*  Y±  jIj-^j-O  for  j>N, 


(5) 


iAx  j=i 


where 


G(u,x)  =  —7—;   u=iAx,  i=m,  m+1,  ..., 

/..2  „2\1/2 


(xz-uz) 
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In  matrix  form,  h=Af.     is  an  upper  triangular  matrix  with  entries 
coefficients  y...      The  calculations  of  h.,  though  they  happen  to  rep- 
resent improper  integrals,  are  straight-forward  and  will  be  described  in 
Chapter  4.   Suffice  here  to  say  that  in  each  summand  of  eq.  (5),  I(x) 
is  expanded  in  a  Taylor  series  around  x=jAx,  retaining  only  the  first 
three  terms  of  the  expansion. 

The  right  hand  side  of  eq.  (3)  is  expanded  in  terms  of  the  discrete 


values  I .  as 

l 


NAx  N 

J     I(x)K(u,x)dx=  /  j  b_.  ,1, 


iAx  j  =  i  1J  J    3 


b..  are  coefficients  to  be  determined.   They  depend  on  the  weighting 
function  W.   Therefore,  from  (5)  and  (6),  \,   which  is  the  N-dimensional 
solution  vector,  is  given  in  matrix  notation  by 

B   is  the  inverse  of  the  matrix  B  whose  elements  are  the  sets  of  the 

coefficents  b..  of  eq .  (6),  and  h  is  an  N-dimensional  vector  with  com- 

ponents  h.,  eq.  (5).   Since  the  matrix  B  is  an  upper  triangular  one,  its 

inversion  can  be  performed  rapidly  by  a  back  substitution,  starting  with 

the  b„  „  element  of  R,  as  follows: 
N,N  u 

VVV'n 

N  (7) 

I.  =  (h.  -/.    b..I.)/b..,  i=N-l,  N-2,  .  .  .  ,  m. 
1   X   j^i+1   ^  J    1X 

The  method  employed  in  this  work  for  solving  eq.  (6)  is  based  on  a 

replacement  of  the  function  I(x)  by  an  interpolating  polynomial  having 
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the  values  I,=I(x,  )  on  a  set  of  points  x,  . 

This  is  a  commonly  used  procedure  employed  in  conjunction  with 
multi-step  methods  which  are  based  on  interpolation.   For  a  multi- 
step  method  of  order  I ,  the  integral  in  eq.  (5)  is  replaced  by  a  repeated 
equally  spaced  quadrature  rule  defined  on  l+l   points.   Thus,  eq.  (3)  is 
written  as  a  system  of  (N-i+1)  simultaneous  equations 

N    (j+l)Ax 
h.  =  y^        J       I(x)K(iAx,x)dx  (8) 

j=i   jAx 


I(x)  is  replaced  by  an  interpolating  polynomial  having  the  values  I,=I(x  ) 
on  a  set  of  points  x  .   In  forward  interpolation  of  multiplicity  I,    the 
interpolating  polynomial  interpolates  I(x)  on  points  x,  ,  x,  1 ,  ..., 
x,   .   In  general,  the  interpolating  polynomial  P  (x)  of  degree  I 
agrees  with  I (x)  and  l+l   suitably  chosen  points,  which  can  be  located  to 
the  right,  to  the  left  or  centered  around  the  given  point  x  ,  representing 
forward,  backward  or  central  interpolations,  respectively.   The  integral 
of  I(x)K(iAx,x)  over  the  interval  [jAx,  (j+l)Ax]  is  then  replaced  by  the 
integral  of  P,  (x)K(iAx,x)  over  the  same  range.   We  will  consider  here 
only  two  cases,  1=1,    and  1=2. 

First,  assume  1=1.      In  the  simplest  case,  which  we  call  the 
"rectangle  rule"  approximation,  I(x)  in  the  r.h.s.  of  eq.  (7)  is  re- 
placed by  I . .   In  this  case, 
J 


(j+DAx 
b. .  =    /    K(iAx,x)dx. 

IT  v 


jAx 
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An  improvement  over  the  rectangle  rule  approximation  for  a  single- 
step  method  (1=1)   would  be,  if  I(x)  in  the  integrals  of  eq .  (8)  are 
replaced  by  the  following  averages  (trapezoid  rule  approximation) : 

I(x)  =  -|(I.+I.+1)  for  jAx<=x$(j+l)Ax.  (9) 

In  this  case,  the  elements  of  the  matrix  g  are: 

m+l)Ax 
b.  .  =  J      K(iAx,x)dx    (j  =  i) 
1,1    iAx 


(10) 


Aj+l)Ax 
b.  .  =  J       K(iAx,x)dx    (j>i) 
1,11    (j-DAx 


b   =0.  (j<i) 

X9  J 


The  question  is:   Is  the  trapezoid  rule  approximation  (eq.  (9))  really 
an  improvement  over  the  rectangle  rule  approximation?   Theoretically, 
the  answer  should  be  yes.   Under  the  rectangle  rule  approximation,  I(x) 
is  assumed  to  be  constant  within  each  strip.   Under  a  trapezoid  rule, 
I(x)  is  assumed  to  be  a  linear  function  of  x  in  each  of  the  integrals 
of  the  r.h.s.  of  eq.  (8),  and  then  replaced  by  its  value  at  the  middle 
points  of  these  integrals.   Nevertheless,  an  application  of  the  trapezoid 
rule  might  lead  to  even-odd  oscillations  in  the  computed  I..   These 
oscillations  are  associated  with  the  presence  of  weak  numerical  instabil- 
ities.  The  weak  instabilities,  in  contradistinction  to  strong  numerical 
instabilities,  can  be  greatly  reduced  by  either  adopting  certain  aver- 
aging procedures  or  by  a  smoothing  method,  and  will  be  discussed  later. 
The  investigations  of  the  numerical  instabilities,  associated  with 

the  employment  of  various  quadrature  rules  can  be  greatly  facilitated  by 
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examination  of  the  elements  of  the  matrix  g   .   We  first  demonstrate  the 
possibility  of  the  presence  of  even-odd  oscillations  in  the  computed 
values  of  I.  under  the  trapezoid  rule  approximation,  and  the  relative 
absence  of  these  oscillations  under  the  rectangle  rule  approximation. 
For  simplicity,  assume  that  W,  the  weighting  function,  is  a  constant, 
equal  to  W=2/tt,  so  that  K(u,x)  =  2x.   Employing  the  rectangle  rule  approxi- 
mation, that  is,  assuming  that  I(x)=I.  in  each  of  the  integrals  of  the 
r.h.s.  of  eq.  (8),  one  obtains  from  eq.  (8), 


b.  .=2j+l  j^i 


b.  .=0  j<i 

i5: 


(11) 


>.  .,  which  are  the  elements  of  the  matrix  R   are  found  to  be 
1,3  u 


i,i   2i+l 
1 


i,i+l     2i+l 


1.  .  =  0  for  -Ui+2 
i,J 


»-li 


Since  J=g  |j,  the  computed  I.  are  seen  to  depend  on  h.  and  h  .  only. 
Thus,  the  rectangle-rule  approximation  will  not  introduce  spurious 
short-range  maxima  and  minima  into  the  computed  I(x)  vs.  x  curves,  if 
the  original  experimental  I(x)  curves  were  reasonably  smooth. 

Next,  consider  the  trapezoid  rule  approximation,  eq .  (9).   Now, 
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h±i±   -  \   (2J+1)         j-i 

b.  .  =  2j  j>i  (12) 

b.  .  =  0  j<i 

i,J 

For  large  i,  the  elements  of  the  matrix  D   are  approximately, 

/. 

for  j=i+2,  i+4,  i+6,  ... 

for  j>i+l,  i+3,  i+5,  ...      (13) 

for  j=i. 


i,j   2i+l 

4 


i,j     2±+l 

3i,j  =  2i+l 

Thus, 

I.  =  h.g.  .  +^7^7  (-h.  .,,+h.  ,,0-h.  .,,+  ...)       (14) 
x    i  1,1   2i+l    i,i+l   i,i+2   i,i+3 

and  the  possibility  of  even-odd  oscillations  in  the  computed  I.  is  self- 
evident. 

Before  getting  into  discussions  of  the  methods  that  might  partially 
eliminate  these  oscillations,  let  us  next  consider  the  case  when  I(x), 
eq.  (8)  is  replaced  by  an  interpolating  polynomial  having  values  on  a 
set  of  points  x,  .   Let  Q (x)  be  the  polynomial  which  approximates  and 
replaces  I(x)  in  eq.  (8). 

This  polynomial  can  be  expressed  in  terms  of  forward,  central  or 

backward  differences.   Assume  1=2.      In  central  difference  interpolation, 

the  interpolating  points  are  x.,  x.  ,  and  x.,n.   In  backward  inter- 

J    3-1      J+1 

polation,  they  are  x. ,  x.    and  x.  ?.   In  forward  interpolation,  they 

J    J —  i      J —  ^ 

are  x.  ,  x.  , ..  and  x.  ,  _  . 
J    j+1      j+2 
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Consider  central  difference  interpolation  (for  i~2 ,  this  method  is 

equivalent  to  the  extended  Simpson  rule  for  numerical  integration) . 

Consider  the  case  of  W=constant=2/ir .   The  interpolating  polynomial  Q(x) 

is  given  by 

(I   -21. +1   ) 

Q(x)  =  I.  +  — :L-i J ±^      (x^jAx)2  (15) 

3  2 (Ax)2 

For  i>>l,  the  elements  of  the  ith  row  of  the  B  matrix  are,  approximately 

b.  .  =  i/3 
1,1 

b.  .  =  4j/3   (j=i+l,  i+3,  i+5,  ...)     (16) 
i » J 

b,  ,  =  2j/3   (j=i+2,  i+4,  i+6,  ...-). 
i>  J 

When  this  matrix  is  inverted,  3.  .,  the  elements  of  the  inverted  matrix, 

i5J 

are  approximately  given  by 

3.  .  =  (-l)1*"1  {(-2+/3)j_i  +  (-2-/3) j_i} 
1,J    i 

Thus,  the  absolute  values  of  3.  .  are  seen  to  grow  exponentially  as  one 

moves  away  from  the  diagonal  elements  of  the  inverted  B  matrix.   As  a 

result,  the  numerical  calculations  are  affected  by  strong  numerical 

instability  which  renders  the  above  method  of  constructing  interpolating 

polynomial  to  be  totally  useless. 

This  strong  instability  is  characterized  as  a  pointwise  phenomenon, 
as  it  involves  error  growth  at  fixed  points  as  the  interval  between  them 
decreases.   In  other  words,  the  multi-step  solution  of  eq .  (8),  when 
I(x)  is  replaced  by  the  above  given  interpolating  polynomial,  diverges 
as  Ax-+0  [13]. 

The  question  is,  how  one  constructs  an  interpolating  polynomial 
Q(x)  that  will  lead  to  numerically  stable  (or  convergent)  method  for 
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calculating  the  solution  vector  I?   This  question  was  answered  in  de- 
tail in  reference  [6].   Suffice  here  to  say  that,  if  the  coefficients 
of  I   in  the  interpolating  polynomial  are  a,  ,  so  that  I(x)  in  r.h.s.  of 
eq.  (8)  is  replaced  by 

a.I.+a.  ,,1..,+.  •  .+a.  ,  'I. ,  . 

i  x  l+l  l+l      x+Z   ±+i 

then  the  moduli  of  all  roots  of  the  characteristic  polynomial 

<Kn)  =  a.n£  +  ai+1**~S«-..^i+£  =  °         <18> 

must  not  exceed  unity,  as  a  condition  for  convergence  of  the  multi-step 
method  of  multiplicity  &1.   In  the  above  example,  in  which  central 
difference  method  was  employed  for  the  construction  of  the  interpolating 
polynomial,  the  characteristic  polynomial  <Kn)  is  given  by 

<Kn)  =  !/3n2  +  V3n  +  1/3 
and  the  root  condition  for  numerical  stability  is  not  satisfied.   In 
reference  [6]  it  was  shown  that,  if  instead  of  central  differences,  one 
employs  forward  differences,  a  convergent  method  can  be  obtained.   (In 
multi-step  method,  employment  of  forward  differences  is  a  necessary,  but 
not  sufficient  condition  for  numerical  stability.)   The  interpolating 
points  for  1=2   are  j,  j+1  and  j+2.   The  root-determining  equation  is 

<K.n)  =  7/3n2  -  2/3n  +  1/3  =  o.  (19) 

Since  the  moduli  of  the  roots  of  eq.  (19)  are  less  than  1,  the  method  is 
numerically  stable. 


■^This  is  a  very  serious  restriction.  As  a  matter  of  fact,  it  is  very 
hard  to  construct  an  interpolating  polynomial  which  will  satisfy  the 
root  condition,  when  £>2. 
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For  i>>l,  the  elements  of  the  ith  row  of  the  B  matrix  (for  constant 
W)  are: 

b.  .  =  7i/3 

1,1 

b.  ,  =  -2j/3        (jfl+1,  i+3,  i+5,  ...) 
i»  J 

b.  .  =  2j/3         (j-i+2,  i+4",  i+6,  ...). 
i  j  J 

We  designate  the  above-described  method  of  constructing  inter- 
polating polynomial  for  I(x)  as  the  Nystrom  algorithm  [12]  for  the 
following  reason: 

In  analogy  to  the  Nystrom  method  for  constructing  interpolating 
polynomial  for  solving  difference  equations  of  first  order,  this  poly- 
nomial is  based  on  forward  interpolation  and  is  applied  to  the  midpoint 
of  each  of  the  summands  in  the  quadrature  rule  approximation  of  the 
integral  equation  (because  it  is  applied  to  the  midpoint,  the  coefficient 
of  VI,  is  zero,  thus  insuring  numerical  stability). 

In  our  computer  program,  this  algorithm  is  employed  in  addition  to 
the  rectangle  rule  and  trapezoid  rule  approximations  for  replacing  I(x), 
eq.  (6)  by  a  set  of  discrete  values  I,  . 

On  the  basis  of  the  above  discussion  on  condition  for  numerical 
stability,  it  becomes  clear  why  the  one-step  methods  are  always  numeri- 
cally stable:   for  the  rectangle  rule  approximation,  the  root-determining 
characteristic  polynomial  (j>(n)  is  of  zero  order.   This  method  is  there- 
fore the  most  stable  numerically.   On  the  other  hand,  the  trapezoid  rule 
method  implies  that  the  characteristic  polynomial  is  of  the  first  order, 
with  root  equal  to  -1.   This  represents  the  border-line  case  for  numeri- 
cal stability. 
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All  these  considerations  apply  strictly  to  the  case  of  W=constant. 
For  a  non-constant  W  we  found  that  the  stability  conditions  are  slightly 
more  relaxed.   Specifically,  for  W  represented  by  a  trapezoid,  it  is 
easy  to  show  that  the  absolute  values  of  the  coefficients  b.  .  decrease 
as  j  increases  as  compared  with  the  case  of  constant  W.   This  fact  is  of 
importance,  particularly  with  conjunction  with  the  trapezoid  rule  method, 
as  it  implies  that  a  non-constant  W  tends  to  eliminate  some  of  the  even- 
odd  oscillations  in  the  computed  I(x),  as  compared  with  the  case  of 
W=constant . 

As  expected,  the  Nystrom  algorithm  leads  to  results  that  are  rela- 
tively free  of  even-odd  oscillations.   This  is  because  3.  .,  the  i,i 

1,1 

element  of  B   matrix,  is  approximately  proportional  to  -r-|  A  P   where 
|a|  is  the  modulus  of  the  largest  root  of  the  characteristic  polynomial 
(f)(n).   Since  X<<1,  the  matrix  elements  3..  are  seen  to  decrease  rapidly 
as  one  moves,  in  a  given  row,  away  from  the  diagonal  elements  of  the 
matrix  B 

As  demonstrated  above,  the  solution  -  vector  \    can  display  undesir- 
able even-odd  oscillations,  particularly  when  |a|,  the  modulus  of  the 
root  of  the  characteristic  polynomial,  is  close  to  one.   This  is  parti- 
cularly the  case  when  trapezoid  rule  is  applied  to  the  solution  of  eq. 
(8) .   These  oscillations  can  be  greatly  reduced  by  smoothing  the  results 
as  developed  in  reference  [6].   In  this  method,  one  computes  first  the 


solution  vector  \=\  from 


-1 


1°  =  B„h  (20) 


B   is  the  matrix  of  b .  .,  calculated  from  eq .  (7,9,  and  10).   The 

"corrected"  solution  vector  \  is  obtained  by  calculating  the  second 
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order  differences  62I .   '=1 .,-21 .   +1 . ,,   for  each  element  I.    of  the 

1      x-1    1     1+1  1 

solution  vector  J    ,  and  adding  y&2I.  to  I .  '  .   For  y ,  a  value  of 

1/6  was  adopted.   In  reference  [8],  y=lM.   However,  we  found,  on  the 
basis  of  numerous  sample  calculations,  that  the  value  y=l/6  is  preferable, 
This  smoothing  procedure  is  repeated,  starting  with  the  solution  vector 
J    ,  to  obtain  a  solution  vector  \         .   Since  the  computations  of  the 
smoothed  solution  vectors  I    and  I    are  very  rapid,  they  are  both 
included  in  the  computer  program. 

In  the  final  print-out,  we  present  solutions  for  the  scattered 
intensity,  based  on  the  rectangle  rule  approximation  and  on  the  trape- 
zoid rule  approximations  (with  and  without  smoothing) .   This  is  done 
for  the  following  practical  purpose:   the  user  of  the  computer  program 
should  recognize  that  some  spurious  small  peaks  in  the  calculated  I(x) 
curves  should  not  be  given  any  physical  interpretation.   Rather,  these 
peaks,  or  oscillations,  in  the  computed  I(x)  curves  result  from  the 
presence  of  weak  numerical  instabilities.   In  order  to  be  sure  whether 
maxima  and  minima  in  the  calculated  I(x)  curves  have  physical  signifi- 
cance, one  should  examine  the  I(x)  curves  obtained  from  the  application 
of  the  rectangle  rule  approximation.   In  the  rectangle  rule  solution, 
weak  numerical  instabilities,  resulting  in  even-odd  oscillations,  are 
usually  absent,  provided  that  the  experimental  I(x)  vs.  x  curves  are 
also  free  of  even-odd  oscillations  (which  might  come  from  noise  effects 
or  from  errors  in  measurements).   Otherwise,  even  with  the  rectangle 
rule  approximation,  these  oscillations  will  be  magnified  in  the  I(x)  vs. 
x  curves.   It  is,  therefore,  very  important  that  the  experimentally 
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determined  I(x)  be  smoothed  before  it  is  read  into  the  computer.   This 
is  true  in  particular  with  respect  to  the  trapezoid  rule  approximation. 

The  solution  based  on  the  Nystrom  algorithm  (multiple  step  method 
of  multiplicity  2,  with  forward  interpolation)  is  more  accurate  than  the 
rectangle-rule  solution,  and  is  relatively  free  of  even-odd  oscilla- 
tions, so  that  corrections  based  on  smoothing  the  results  need  not  be 
applied  to  this  solution.   The  solution  based  on  this  method  should  be 
employed  for  display  of  desmeared  intensities. 

The  various  algorithms  and  their  effect  on  the  calculated  I(x)  are 
discussed  in  Chapter  6,  where  sample  calculations  employing  various 
algorithms  are  presented  and  discussed. 
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CHAPTER  4 

1.   DESCRIPTION  OF  CALCULATIONS 

A.   Preparation  of  Data  which  Characterize  the  Slit 
Weighting  Function 

The  following  slit  weighting  functions  are  employed  in  this  pro- 
gram:  1)   Weighting  function  is  a  constant;  2)   Gaussian  Function;  3) 
Data  Function;  and  4)   Trapezoid  Function. 

For  Gaussian  distribution,  W(r)  is  computed  at  equal  intervals 
equal  to  A.   Thus,  the  program  computes  a  set  of  discrete  values: 

W.  =  exp(-i2A2/G) 
G  is  read  from  the  input  data  card,  and  A  is  arbitrarily  chosen  to  equal 
0.1/G.   The  subscript  i  assumes  values  of  0,  1,  ...,  20. 

W(r)  is  represented  by  data  function.   A  set  of  cards  giving  the 
various  pairs  of  W(r)  and  r  are  read  into  the  computer.   The  data  are 
represented  as  ri(=0),  r£,  r3,  ...,  r  ;  W(ri),  W(r2) ,  W(r3) ,  ...,  W(r  ). 
W(r.)  are  given  at  equally  spaced  intervals.   The  first  pair  of  data 
are  stored  as  ?2   and  W(r2) ,  respectively.   r^  is  set  as  equal  to  zero, 
and  W(r^)  is  computed  according  to 


W(0)=W(r1)=W(r2)+[W(r2)-W(r3) ]r/(r3-r2) . 

There  are  m  pairs  of  data.   If  m<N,  the  weighting  function  W(r)  for 

r>r   is  given  a  constant  value  which  is  equal  to  W(r  ).   Thus,  W(r,  ) 
m  n  m        '    k 

W(r  )  for  m<k<N. 
m 


W(r)  is  represented  by  a  trapezoid,  as  follows 
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W(r) 


W(r)=l     r<S. 


w(r)"s^i  Sl*rSSz 


W(r)=0 
W(-r)=W(r) 


r>S, 


S.   Evaluation  of  the  Integral 

oo 

/%,    v           NAx  y,    v 
xl(x) dx  _  /•   xl(x) 
(x2-u2)1/2    "#.A  [x2-(iAx)2]1/2 
u               ,/iAx 


(21) 


since  l(x)=0  for  x>NAx.   The  integral  h(u)  =  h(iAx)  is  computed  for  each 

value  of  i  by  the  subroutine  CGEN  which  is  entered  with  a  given  value  of 

i.   To  evaluate  eq.  (21)  we  subdivide  the  integral  into  N-i-1  integrals 

as  follows: 

N      (j-i)Ax 


i-E  ../- 


h.  = 


xl(x) 


jAx   [x2-(iAx)2]1/2 


dx 


(22) 


I(x)  is  expanded  in  a  Taylor  series  around  x  =  jAx 
N    (j+l)Ax 


^S/  iv 


j=i  jAx 


I., -.-I.    , 
[x-(j+l)Ax]   3^  J  +  j   [x-(i+l)Ax]2 


(I   ..-2I.+I.  .) 

3+1   3   3-1  >      x 

(Ax)2       [x2-(iAx)2]1/2 


(23) 


dx 
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N 
h.  =X^y..I.;    I.=0  for  j>N.  (24) 


%    ^ 


In  eq.  (24)  and  in  following  equations,  I.  =  I(jAx) ,  etc.   y.  .  form  an 

3  i  >3 

upper  triangular  matrix  Q.   The  matrix  elements  y.  .  have  different 

i  »3 

representations  depending  on  the  following  three  cases:   (1)   j-i=0,  (2) 

j-i  even  and  (3)   j-i  odd.   In  addition,  in  order  to  close  the  system 

one  has  to  consider,  separately,  the  cases  of  j=N  and  j=N-l. 

Sometimes  the  high  accuracy  in  calculating  h.  is  not  required.   One 

can  simplify  calculations  and  save  computer  time  by  retaining,  in  the 

r.h.s.  of  eq.  (23),  only  the  first  term  in  the  Taylor  expansion  of 

I(x).   This  option  is  controlled  by  the  parameter  MM,  read  with  the  data 

cards.   Since  in  this  case  the  results  are  simple,  we  spell  them  here. 

They  might  be  of  use  for  a  fast  check  of  the  results: 

y.  .=1/2  {[(j+l)2-i2]1/2  -  [(j-1)2  -  i2)1/2}Ax  for  N-l>j>i 
1  >  J 

y.  .=1/2  (2i+l)1/2Ax  for  j=i 

1  5  3 


y.  .=0  for  i<i. 
i>3 


i!/2  _  r^_^2_,2ll/2- 


For  j=N,  y.   =1/2  { [ (j2-i2) ]x/   -  [ (j-1) 2-i2]X/ Z}Ax. 
1 » J 

C.   Evaluation  of  the  Kernel  of  the  Integral  Equation 
for  Various  Weighting  Functions 

1.   Preliminary  Consideration 

We  wish  to  determine  the  coefficients  b.  .in  the  set  of  simul- 

i>3 

taneous  equations 


JN 

E 


3  =  1 


b.  .1.  =  h.  (25) 

i,3  3     i 
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given  that 

NAx 


h.  =  J     I(x)K(iAx,x)dx  (26) 

iAx 


/-[x2-(iAx)2]1/2 
K(iAx,x)  =  2x  /  W(t) dt       (27) 

[x2-(iAx)2-t2]1/2 


The  integral  of  eq.  (27)  is  replaced  by  the  sum  of  integrals 


j=i  jAx 


(j+DAx 

I(x)K(iAx,x)dx  (28) 


We  assume  that  in  eq.  (28),  I (x)  can  be  replaced  by  its  average 

value  within  each  range  and  taken  outside  the  integral  sign 

provided  that  I(x)  is  continuous  on  [jAx,  (j+l)Ax]. 
Define 

[(j+l)2-i2]1/2Ax   t 

r         C  C       W(r)     dr 

(jz-iz)  'Ax  o 


b.  .  are  identical  to  b .  .  only  when  I(x),  eq.  (28),  is  equal  to  I.. 
Otherwise,  b.  .  will  be  computed  from  b.  .,  according  to  the  algorithm 
employed  in  replacing  I(x),  eq.  (28),  by  a  set  of  discrete  values  I,  . 


Introduce 

Rn  =  (jz-i2)"/"Ax,   R0  =  [(j+l)z-iz]^Ax. 


2_42n1/2Av    d   ^  r^a-^2_,2i1/2, 


From  eq.  (29),  after  changing  order  of  integration,  we  obtain 


35 


R2         R2  R1  R2  (30) 


.        =    f     drW(r)  f     t   ;  ,0   dt  +  f       drW(r)  f       — 


1/2 


dt 


Integrate  over  t,  to  obtain 

R2  Rl 

b.  .  =  f      (R22-r2)1/2W(r)dr  -  j     (R12-r2)1/2W(r)dr .   (31) 

o  o 


Write 

b..  =  B.  ,,,-B.  ..  (32) 

To  compute  B.  .  the  FORTRAN  program  calls  a  Function  subprogram 

B(I,J)  for  each  pair  of  values  i  and  j.   This  procedure  is  employed  in 

order  to  prevent  storing  large  matrices.   In  another  program  which  is  to 

be  employed  when  N<140,  the  elements  b.  .  are  stored  in  a  matrix 

i,J 

array,  and  the  FORTRAN  program  calls  a  subroutine  B(I,J).   The  function 
B(I,J)  or  the  subroutine  B(I,J)  calls  other  subroutines  for  evaluation 
of  the  integrals  in  eq.  (31).   The  calculations  for  each  pair  of  vari- 
ables i  and  j  are  performed  twice.   The  subroutine  that  computes  eq. 
(31)  is  entered  with  a  variable  R=R2  the  first  time,  and  with  a  variable 
R=Rl  the  second  time. 

Both  the  Gaussian  and  the  Data  function  representation  for  W(r)  are 
treated  in  the  same  way,  except  that,  for  the  Gaussian  W(r) ,  W(r)  are 
first  computed  by  the  Main  Program.   The  computation  of  B.    _  and  of 
B.  .  for  the  Gaussian  and  Data  function  representations  of  W(r) 
are  performed  in  the  subroutine  WGHT23.   For  the  trapezoid  W(r) ,  these 
computations  are  performed  in  the  subroutine  WGHT4.   For  W(r)  =  constant, 


b±(.  -Jsj+i). 
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2.   Computation  of  eg.  (31) 

a)   W  is  given  by  a  set  of  data.   Subroutine  WGHT23  is  entered  the  first 

1  /  7  1  /  9 

time  with  R  =  [(j+l)2-i2]   Ax,  the  second  time  with  R  =  (j2-i2)    Ax. 

Since  R  is  generally  a  non- integer,  we  locate  an  integer  R,  such  that 

Calculate  W(R)  by  linear  interpolation  as  f ollows : 

W(R    )-W(R  ) 
W(R)  =  W(R  )  +    ^     —  (R-R  )    (r  >R) 

or 

W(r  )-W(r    ) 

W(R)  =  W(r  )  '.+  - ^-±-   (R-r  )    (r  $R)  .      (34) 

m      r  -  r   -       m      m 

m    m-1 

W(r  )  and  r   is  the  last  pair  of  data  read  into  the  computer.   Each  of 
mm 

the  integrals,  eq.  (31),  is  represented  by  a  sum  of  smaller  integrals. 

B .  . , .  is  written  as 
i,l+l 

ft    rk  R2 

B±  +1  =  /    /   (R?2-R2)1/2  W(r)dr  +  /  (R02-r2)1/2W(r)dr .    (35) 

k=1   Vl  r£ 


L2 


In  each  of  the  integrals  in  the  summation,  eq.  (35),  W(r)  is  approxi- 
mated by  a  linear  dependence  on  r  within  each  strip  (r,   ,r,).   The 
individual  integrals  are  expressed  as  follows: 

x 

k       f-n    2       2\1/2tt/'    \j  k      k-l    /  r-n    2       2n1/2j 

(R0z-r^)        W(r)dr  =  — ^ /  (R0    -rz)         dr  + 


/ 

-1 

(36) 


2  '        i^'."«  2  J  v    2 

rk-l 
rk-l 


rk 


j  -w         r 

k     k-l        / 
:k~rk-l       J 


W,  -W,     ,  /  ,  r  +r 

(R  2-r2)1/2    (r-  -*-£%    dr, 


k-l 
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In  eq.  (36),  W  =W(r.  ) ,  etc.   In  eq.  (36)  the  last  term  is  an  improvement 

rC      K. 

over  the  trapezoid  rule  approximation. 

Performing  the  integrations  and  rearranging  terms,  one  obtains 


'i>J+1  >  I    £  M»t.1-»k«k+i)  {V  sln"1  (f)    +  Wh2-^2^2 


+  354  [r22-(^)2i3/2|  +  jVi-  7^7  <VV-M 

R„-£A  2 


X  |  R^sin"1    (M)    +  £A[R22-(£A)2]1/2  +   3^   [R^   -    (£A)2]3/2 
ttR   2  R  +U 

+  ._jWr2  +  w^____  (wR2-w£)}. 


In  eq.  (37),  r  =kA,  A=r,-r,   .   In  deriving  eq.  (37),  use  was  made  of 
the  fact  that  r  =0,  by  definition.   In  deriving  the  last  two  terms  of  eq, 
(37),  (£-l)W  A  and  R„W  A  were  approximated  by  £W  A .   Similar  results  are 

/V  i—  Xj  Xj 

obtained  for  R..  replacing  R„ .   As  a  fast  check,  assume  a  constant  W(r)  = 

W.   In   this  case,  b.  .  =  -t-(R  2-R  2)  ,  which  also  obtains  directly  from 

1,3   4  2   1 

eq.  (31). 

b)   Trapezoid  Weighting  Function.   Let 

W(t)  =  1  (t<s1) 

W(t)  =  (S2-t)/(S2-S1)     (S1<t<S2)  (38) 

W(t)  =  0  (t>s2) 
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Substitute  (39)  into  (30)  to  obtain 


1)   R2>S2 


R2_s  2)1/2  +    2s.n 


'i,j+l   "lv^2    1 


-S1(R22-S12)1/2+R2 


l> 


•  -if- 

2   "in   \R )'    "2V"2 


S0(R  2  -  S  2)1/2 


".  -l/S2\    .  -l/Ml    S2 
sxn  y—xn   ^Js^ 


(39) 


+|  [(r2-S  2)3/2  _  (R  2_s  2)3/2] 


3  LV  2    2 


2    1 


S2"S1 


2)   S1<R2<S2 


1/2 


i.1+1  =  W  "  Sl  >    +  V  S±n  \TJ  +     2   V  "  Si(R2Z-Sl  > 


'©•  i 


2  c  2Nl/2 


-  R,/  sin 


.  -1  /Sl\   S2 


2  (R22-S;L2)3/2 

3  (S2-S1) 


(40) 


3)   R2<S1 


Ei,j+1  =  2  R2, 


(41) 


Similar  expressions  are  obtained  for  B.  .,  with  R,  replacing  R„ . 

D.   Calculations  of  b.  . 

i,3 

For  the  rectangle  rule  approximation,  b.  .  =  B.  ..   For  the  trape- 

i,3    i,3 

zoid  rule  approximation,  we  calculate  first  solution  vector  \  for 

the  rectangle  rule  approximation, 

I (o)  -  B_ah 


where  R   is  the  inverse  matrix  of  the  elements  b.  ..   Then, 

i,3 
t (trapezoid)  _  n_lT (°) 

the  elements  of  the  £)  matrix  are 


(42) 


(43) 


D.  .  =  D.  ...  =  1/2 
1 ,  X     X  ,  x+1 

D.  .  =  0  j<i+l  and  i>i. 
i,3 
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Therefore,  B  =  D  B-   F°r  the  Nystrom  algorithm  (I(x)  extrapolated  forward), 

we  employ  the  following  approximate  transformation:   Form  matrix  f"  of 

elements  f.  .as  follows: 
i,3 

£.  ,-Ib.  . 

1,1   3  1,1 

f,  ,   =   f  b    (j=i+l,  i+3,  ...)  (44) 

-1-  j  J    ->   1»3 

fl,J-fEi.3  ^+2-1+*. ■■■■•)- 

Equations  (45)  are  based  on  the  following  arguments:   To  insure  numeri- 
cal stability,  the  point  of  reference  for  interpolating  polynomial  must 
be  located  in  the  middle  of  the  integrand.   Therefore,  eq .  (28)  is 
replaced  by  N  (j+l)Ax 


u  =Z^  J    I(x)l 


h_.  =  /  j  J    I(x)K(iAx,x)dx  (45) 

j=i,i+2,  ...  (j-l)Ax 

Q(x) ,  the  interpolating  polynomial  which  replaces  I(x),  is  given  by 

Q(x)  =  I  +  (I  -21  +1+Ij+2)  (x-jAx)2/2(Ax)2.        (46) 

Therefore,  the  calculations  of  I.  must  start  with  i=m+l,  since  with  this 

l 

algorithm,  calculation  of  I  will  always  be  numerically  unstable,  as 
then  the  root  condition  (eq.  18)  will  not  be  satisfied.   To  close  the 


system, 


and  of  course, 


f.  ,  =  2/3  b.       (j=N,  j-i  even) 

f.  .  =  1/3  b.  .     (j=N,  j-i  odd)  (47) 

i,J        i,3 

f.  .=2/3  b.         (j=N-l,  j-i  even) 
i,3      1,3 


b .  .  =  f .  . ,  B  =  R 

1,3     i,3 
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Therefore, 

T(Nystrom)  ^  1--I1 

i  =  h   n- 

E.   Solution  Checking 
In  order  to  be  able  to  make  a  direct  check  on  the  accuracy  of  the 
original  calculations,  we  recalculated  I(x)  by  numerical  integration, 
from 

NAx  (N2-i2)1/2Ax 

cy   o  f   yi(y)w(y2  -  (iAx)2)1/2     if  w(t)  (i[t2-(iAx)2]1/2)dt. 

i  ~  J a     r  2    ,.A  ^71/2     dy  =  •/  (48) 

lAx    [yz  -  (iAx)z]  o 


.2_,,a^2m1^  .. 


I( [t  -(iAx)z) ]    )  is  calculated  from  the  computed  I  by  quadratic 
interpolation. 
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CHAPTER    5 
FORTRAN    PROGRAM 
C  CALCULATION    FOR     SLIT    CORRECTION     IN     SMALL    ANGLE    X-RAY    SCATTERING 

C  PROGRAM    AVOIDS    MATRIX    STORAGE 

DIMENSION    XI  (4  00  ),  YH400  ),XX(4  00  ),Z  1(4  00  ),YK(  40  0  ),YL(  40  0),YY(  4  00), 
1  YM(  400), U(  400) 

DIMENSION    LW(  10  ),WM2(  10  ) 

COMM0N/BLOCK1  /A,  N,  I  START,  NN,  MD,  MM 

C0MMON/BI  0CK2/RN(  400  ),W(  400  ),V(  400  ),C(  400  ),DL,G 

COMMON    /BL0CK3/    KB,KL 

COMMON    /BL0CK5/S1,  S2 

DATA    PI/3.1415926/ 

DATA    IPR,IRD/6,5/ 
C  READ    DIRECTIVE     I NSTPUCTI ONS    AND    SET    NUMBER    OF    COMPONENTS     IN     INPUT    VECTOR 

WRITE(  IPR,1  ) 
C  AT    SOME     POINT    THE     INTERVAL     BETWEEN     INPUT    DATA     IS    GREATER    THAN    AT 

C  INITIAL    ANGLES. THAT    POINT    IS    GIVEN    BY     * IENTR'. THE    NUMBER    OF    SUBDI- 

C  VISIONS    BETWEEN    CONSECUTIVE     INPUT    DATA    FROM    THIS    POINT    IS    GIVEN    BY 

C  'NSUB1 . 

C  N    IS    THE    NUMBER    OF     COMPONENTS    OF    THE     INPUT    VECTOR    AFTFR     INTERPOLA- 

C  TION.NN    IS     THE    NUMBER     OF     INPUT    DATA    READ     INTO    THE    MACHINE. 

C  I  START    IS    THE    RATIO    OF     THE    LOWEST    ANGLE    AT    WHICH    MEASUREMENT     IS 

C  PERFORMED    TO    THE     INTERVAL     BETWEEN    TWO    CONSECUTIVE    ENTRANCES. IT     IS 

C  GIVEN    AS    INTEGER 

READ(  IRD,  1  0  )    NN,  ISTART,  IENTR,  NSUB 

SSUB-NSUB 
1 0  FORMAT    (4110  ) 

1  F0RMAT(  1H1     //    10X,15H    SPECIFICATIONS) 

N-NSUB*(  NN-IENTR*I START )*NN* ISTART 
WRITE(  IPR,30)     N,NN,  ISTART 

30  FORMAT(  1X.2  3HN,  TOTAL    RANGE    OF    DATA-     , 13, 4X, 24HNN, NUMBER    OF    DATA     RE 
IAD-     , 13, 4X.24HI START, STARTING    POINT,"     13) 

IF(  IENTR.  EQ.O  )    GO    TO    36 
WRITE(  IPR,35  )     IENTR,  NSUB 

35  F0RMAT(  1 X, 33HENTRANCE    POINT    FOR     SUBDIVISION,"     , 13, 4X, 27HNUMBER    OF 
1SUBTABULATI0NS,  »     ,13) 

36  K0"0 

31  K0"KO*1 

C      DATA  FOR  WIDTH  CORRECTION. WM2  IS  THE  SECOND  MOMENT  AND  LW  IS  THE 
C      FIRST  POINT  AT  WHICH  THE  WIDTH  CORRECTION  APPLIES.  SEVERAL  WIDTH 
C      CORRECTIONS  CAN  BE  PERFORMED  IN  SUCCESSION. 
READ(IRD,32)  LW(  KO  )  ,  WM2(  KO  ) 

32  FORMAT(  I10.F10.5  ) 

IF(  LW(  KOKEQ.O  )    GO    TO    33 
WRITE(  IPR,37  )    LW(  K0),WM2(  K0) 
3^  FORMAT(  5X,3  0HWIDTH    CORRECTION    APPLIED    FROM    ,  15, 8X, 1 5HSEC0ND    MOMENT 

1  ■     ,F10.5) 
GO    TO    31 

33  KO-KO-1 
WRITE(  IPR,200) 

3  8  FORMAT(/) 

I0-ISTART*1 

M"N-1 

N1"N*1 
C        SELECT  THE  WEIGHTING  FUNCTION 

READ(IRD.IO)    MD,IDC, INV,MM 
GO    TO    (40,60,100,230  ),MD 
C  WEIGHTING    FUNCTION    IS    CONSTANT 

40  WRITE(  I  PR, 50  ) 

50  F0RMAT(  20 X,31 HWEIGHTING    FUNCTION    IS    CONSTANT    // ) 
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c 
c 
c 

60 
70 


80 
90 


C 

100 

110 

120 


GO  TO  2  60 
PREPARATIONS 


FOR  GAUSSIAN  WEIGHTING 


130 
140 


C 

150 


160 

C 

170 

180 


190 


200 
21  0 

220 

C 

230 

24-0 


REPRESENTATION 
EQUAL  TO  0.02 


OF 


FUNCTION 
WEIGHTING 


FUNCTION     ARE    TERMI' 


,  F7.  3,  4H    G-     ,F7.1//) 


OF    DATA    FOR     THE    DATA    WEIGHTING    FUNCTION 


DATA    FOR    GAUSSIAN 
NATED    WHEN    W(  I  )     IS 
READ(  IRD,70)    G 
FORMAT(  2F10.2  ) 
DL«SORT(  0.01  *G  ) 
DO    80    I  -1  ,21 
EI-I 

RN(  I  )-(  EI  -1  .0  )»DL 
W(  I  )  -EXP(  -RN(  I  )*#2/G  ) 
WRITE(  IPR.90  )    DL,G 

FORMATC  4X,19HW    IS    GAUSSIAN,     DL  = 
KBM 
KL»21 
GO    TO    2  60 

PREPARATIONS 
KM 
K»K*1 

READ(  IRD,120  )    RN(K),W(K) 
FORMAT    (F9.3,F11.5) 

IF    (  RN(  K).GE. 0.00000001  )    GO    TO    110 
IF    (  W(  K).GE. 0.00000001  )    GO    TO     110 
DL-RN(K-1   )-RN(K-2) 
EN-N 
IF    (  RN(  2).E0.0.  )    GO    TO    130 

FIND    W(l  )    BY    EXTRAPOLATION 
KB-1 

RN(  1  )  -  0  . 

W(  1  )-W(  2)*(  W(  2  )-W(  3)  )*RN(2  )/(  RN(  3  )-RN(  2  )  ) 
GO    TO    140 
KB«2 

IF    (  RN(  K-l  ).LE.EN)    GO    TO    150 
KL«K-1 
G6    TO    1  70 

EXTEND    THE    WEIGHTING    FUNCTION    TO    THE    END    OF    THE    RANGE    OF    COMPUTATIONS 
KL-EN/DLM  .5 
DO    160     I"K,KL 
RN(  I  )-RN(  1-1  )*DL 
W(  I  )«W(  K-l  ) 

PRINT    OUT    THE    WEIGHTING    FUNCTION 
WRITE( IPR,180) 

FORMAT    (50X,23H    W    IS    A    DATA    FUNCTION.     //34X,13H    GIVEN    DATA.  ,//41  X, 
2     2H    K,19X,2H    R,19X,2H    W// ) 
DO    210    K-KB.KL 

WRITEC  IPR.190  )    K,RN(K),W(K) 
FORMAT    (40X,I3,14X,F9.3, 12X,F1  1.  5/  ) 
SK-K/5 
XK-K 

IF    (  XK/5.0.NE. SK )    GO    TO    210 
WRITE(  IPR.200) 
FORMAT    (/) 
CONTINUE 
WRITE(  I  PR, 22  0  ) 
FORMAT    (  1  HI  ) 
GO    TO    2  60 

PREPARATIONS    FOR    THE    TRAPEZOID    WEIGHTING    FUNCTION 
READ(IRD,70)    SI  ,  S2 
WRITEC  IPR,240)     SI  ,  S2 
FORMATC  4X.19HW    IS    TRAPEZOID,  SI  -     ,F7.3,5H    S2-     ,F7.3//) 
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DL-1.0 

DO    250     IM,N1 
250  RN(  I  )-I  -1 

C  READ    IN    INPUT    VECT6R 

C  IF    BLANK, STOP    THE    CALCULATIONS 

C  READ    ONLY    N    COMPONENTS 

260  READ(  IRD.270)    (  XI(  I  ),  I  -  1  ,  8  ) 

270  F6RMAT  (8F10.0) 

IF  (  XI(  1  ).LE. 0.00000001  )  GO  TO  800 

READ( IRD,270 )  (  XI(  I  ),  I -9 , NN  ) 
C      SHIFT  THE  INPUT  DATA  BY  I  START 

D6  265  II *1,NN 

T"NN*1  -II 

II -1*1  START 
265  XI(  ID-XK  I  ) 

C  DATA     INTERPOLATION    BY    QUADRATIC    FORWARD     INTERP0LATION . IF     IENTR    IS 

C  EQUAL    TO    ZERO, DO    NOT    INTERPOLATE 

IF( IENTR)    276,276,271 

271  S-l  ,0/(  1  .0*SSUB) 
N0»NN*ISTART-2 
N2«NN*I START 
JENTR-IENTR 
K-JENTR 

DO    273    I»JENTR,N0 
XX(  K)»XI(  I  ) 
DO    272    J»1,NSUB 
T-J 
K«K*1 
27  2  XX(  K)-XI(  I  )»T*S#(  XI (  1*1  )-XI(  I  )  )*(XI(  I*2)-2.*XI(  1*1  )*XI(  I)  )*T#S 

1  *(  T*S-1  .  )/2.0 

273  K-K*l 

XX(  K)»XI(N2-1  ) 
DO    274    J*1,NSUB 
T»J 
K*K*1 

274  XX(  K)-XI(  N2-1  )*T*S*(  XI(  N2  )-XI(  N2-1  )  ) 
XX(  K*l  )-XI(  N2) 

DO    275    I-JENTR,N 

275  XI(  I  )»XX(  I  ) 

276  CONTINUE 

C      CORRECT  INPUT  DATA  FOR  WIDTH, IF  NEEDED 

IF(KO.GE.l)  CALL  WIDTHC  XI  ,  LW,  WM2,  KO  ) 

DO  290  I«IO,N 

CALL  CGEN  (  I  ) 
C        MULTIPLY  MATRIX  C  BY  THE  INPUT  VECTOR 
C        CALCULATE  AND  STORE  AS  V( I  )  THE  ITH  ROW  OF  MATRIX  C 

C(  I  )-0.0 

DO  290  J-I,N 
290    C(  I  )-C(  I  )*V(  J)*XI(  J  ) 

C        MATRIX  INVERSION  BY  BACK  SUBSTITUTION,  ZEROTH  ORDER  APPR0XIMATI9N 
31 0    DO  320  I-IO.N 
320    YY(  I  )-C(  I  ) 

A-l  .0 

DO  340  Kl  -10, N 

K-I0*N-K1 

DO  340  I"I0,X 

IF  (I.EO.K)  GO  TO  330 

YY(  I  )-YY(  I  )-YY(  K)*B(  I  ,  K  )/B(  K,  K  ) 

GO  TO  340 
33  0    YY(  I  )-YY(  I  )/B(  1,1  ) 


44 


340  CONTINUE 

C  OPTIONAL     PRINT-OUT 

IF    (IDC.E0.0)    GO    TO    390 

IF(MM.NE.O)  WRITE(  IPR,380  ) 

WRITE(  I  PR, 350) 
350    F0RMAT(  ///8X,  1  1  HC-VALUES. .  /// ) 

DO  360  I«IO,N 

II-I-l 
360  WRITEC  IPR.370)     1 1 ,  C(  I  ) 

370  FORMATC  15X,I4,E15.5  ) 

380  FORMATC  8X,35HL0W    6RDER     APPROXIMATION     IN    C- VALUES/ ) 

C  MATRIX    INVERSION, TRAPEZOID    RULE    APPROXIMATION 

390  DO    400     I"IO,N 

400  YI(  I  )*YY(  I  ) 

YI<  N)«YY<  N)/2.0 

D0    420    Kl "10, N 

K-I0*N-K1 

D0    420    I-IO,K 

IF    (I.EO.K)    G0    T0    410 

YI(  I  )»YI(  I  )-YI(  K)*D(  I,K  )/D(K,K  ) 

G0    T0    420 
410  YI(  I  )-YI(  I  )/D(  1,1  ) 

42  0  CONTINUE 

C  MATRIX    INVERSI0N,NYSTR0M    ALGORITHM    (FORWARD    INTERPOLATION) 

A-0. 
DO    430     I»IO,N 

43  0  YM(  I  )-C(  I  ) 

YM(  N)»YM(  N)/2.0 
DO  450  Kl  "10,  N 
K«I0*N-K1 
DO  450  I»I0,K 
IF  (  I.EO.K)  G0  TO  440 
YM(  I  )-YM(  I  )-YM(K)*F(  I,K)/F(K,K) 
GO  TO  450 
440    YM(  I  )"YM<  I  )/F(  I,  I  ) 

450  CONTINUE 

C        SHIFT  OUTPUT  VECTOR  YM  BY  ONE. SET  YM(  1  )  TO  ZERO, SINCE  IT  CANNOT 
C      BE  COMPUTED  WITH  THIS  ALGORITHM. 

DO  451  II -10,  N 

I-I6*N-I1 

451  YM(  1*1  )-YM(  I  ) 
YM(  I0)»O. 

C        APPLY  LOCAL  PERTURBATIONS  TO  REDUCE  EVEN-ODD  OSCILLATIONS 

470    DO  480  1*2, M 

480    YK(  I  )-(  YI(  1-1  )-2.0*YI(  I  )*YI(  1*1  )  )/6.0 

YK(  1  )-0.0 

DO  490  I-l.M 
490    YK(  I  )-YI(  I  )»YK(  I  ) 

DO  500  I»2,M 
500    YL(  I  )-(  YK(  1-1  )-2.0*YK(  I  )*YK(  1*1  )  )/6.  0 

YL(  1  )-0.0 

DO  510  1-1  ,M 
51  O    YL(  I  )«YK(  I  )*YL(  I  ) 
C        PRINT  OUT  RESULTS 

54  0    WRITE(  I  PR, 545 ) 

545  FORMAT(  //    8X,12HINPUT    VECTOR, 2X,     14HRECTANGLE    RULE, 2X, 9HTRAPEZ0 ID, 

1 2X, 1 5HTRAPEZ. SMOOTHED,  2X,  1 4HYK(  I  )    SMOOTHED,  4X,  2 1HNYSTR0  M    INTERPOLA 
2TI0N) 
WRITE(  I  PR, 550) 

55  0  FORMAT(  /3X,  1  HI  ,  6X,  5HXK  I  ),  10X,  5HYY(  I  ),  9X,  5HYK  I  ),  10X,5HYK(  I  ),10X, 
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15HYL(I  ),1  0X,5HYM(  I  )/ ) 

J-0 

DO    560    I'ie,N 

II-I-l 

J*J*1 

IF    (J.LE.10)    GO    TO    560 

J-l 

WRITE(  I  PR, 38  ) 
560  WRITE<  IPR.580)     1 1  ,  XI  (  I  ),  YY<  I  ),  YI(  I  ),  YK<  I  ),  YL(  I  ),  YM(  I  ) 

580  FORMAT    (I4,6E15.6) 

C  IF     NO    BACK    CHECK,  FIND    OUT    WHETHER     THERE    ARE    MARE    DATA    TO    BE    READ 

C  IF    NOT, GO    TO    STOP 

IF    (INV.EO.O)    GO    TO    260 
C  PREPARATIONS    FOR    THE    BACK    CHECK 

C  THE    OUTPUT    VECTOR    TO    BE    CHECKED    IS    PUT    INTO    YI    ARRAY. 

C  INVERSE    CALCULATIONS    BY    EXTENDED    SIMPSON    QUADRATURE 

IF    (MD.EO.l)    GO    TO    7  30 

IF    (MD.EQ.4)    GO    TO    760 
C  DATA    PREPARATIONS    FOR    THE    DATA    WEIGHTING    FUNCTION    FOP     INVERSE    CHECK 

EN-N 

KL-EN/DL*1 .5 

DO     610     1-2, N 

S-I 

DO    600    K-2.KL 

1 1  -  RN(  K  ) 

IF    (  I.LE.II  )    GO    TO    610 
600  CONTINUE 

61  0  U(  I  )»(  W(K)-W(K-1  )  )*(  S-RNCK-1  )  )/DL*W(K-l  ) 

U(  1  )  -W(  1  ) 
620  DO    630    1*1 ,N1 

63  0  RN(  I  )-I  -1 

C  DETERMINE    OUTPUT    VECTOR    ON    WHICH    BACK-CHECKS    ARE    PERFORMED 

DO    590    I-1,N 
C  DETERMINE    OUTPUT    VECTOR    ON    WHICH    BACK    CHECK    IS    PERFORMED. YI( I ) 

C  OF    THE    RIGHT    HAND    SIDE    OF    NEXT    CARD    CAN    BE    REPLACED    BY    YY(  I ) , 

C  YK(  I  )  ,  YL(  I  )  ,  OR    YM(  I  )  . 

590  YI<  I )-YI(  I  ) 

DO    680    I-Ifl,N 

S-I 

ZI(  I  )-W(l  )»YI(  I  )/l  .5 

DO    680    J-2,N1 

T-J 

R-SORT(  (S-l  .0  )**2*(  T-1.0  )**2  ) 

Rl -Nl *1 

IF    (R.GE.R1)    GO    TO    680 

DO    640    K-I0.N1 

IF     (R.LT.RN(K))    GO    TO    650 
640  CONTINUE 

Yl  -YK  N)*(  YI(N)-YI(  N-l  )  )*(R-RN(N  )  )/(  RN(  N  )- RN(  N-  1  )  ) 

GO    TO    660 
65  0  Yl-YK  K-l  )♦(  YI(  K)-YI(K-1  )  )*(R-RN(K-1  )  )/(  RN(  K  )-RN(K-l  )  )♦(  YI(  K*l  )-2. 

20*YI(  K)*YI(  K-l  )  )*0.5*(  R-PN(K-1   )  )*(  (  R-RN(  K- 1  )  )/(  RN(  K)-RN(K-1  )  )-1.0) 
3/(  RN(K)-RN(K-1  )) 
660  IF    (  MOD(  J, 2 ).EQ.0 )    GO    TO    670 

ZI(  I  )-ZI(  I  )*U(  J  )*Y1*4.  0/3.0 

GO    TO    6  80 
670  ZI(  I  )-ZI(  I  )*U(  J  )»Y1#8. 0/3.0 

68  0  CONTINUE 

C  PRINT    OUT    THE    ORIGINAL    AND    THE    RECONSTRUCTED    OUTPUT    VECTORS 

WRITE(  IPR,690) 
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690 
700 
710 


720 
72  5 


72  7 
726 
C 

730 
74  0 


C 

760 


77  0 

780 
790 


C 

800 
81  0 


F0RMAT<  1H1    //    20X.11HBACK    CHECKS    /20X, 5HINPUT, 5X, 1 OHREC OMPUTED ) 

WRITE<IPR,710  ) 

FORMAT(  //15X.27HI  XI(I)  ZI(I)/) 

J-0 

DO    720    I-IO.N 

II-I-l 

J-J*l 

IF    (J,LE.10)    GO    TO    720 

J-l 

WRITE(  I  PR, 38) 

WRITE(  IPR.725  )     1 1 , XI (  I  ),  ZI( I  ) 

FORMATC 12X,I5,2E15.5 ) 

IF(MD.EQ.4)    G6    TO    726 

DO    727    1-1,50 

EI-I 

RN(  I  )-(  EI  -1.0  )*DL 

GO    TO    26  0 

DATA    PREPARATI6NS    F6R    CONSTANT   WEIGHTING    FUNCTION    FOR     INVERSE    CHECK 
DO    740    1*1, Nl 
U(  I  )-l  .0 
W(  1  )  -1  .  0 

go  to  620 
data  preparations  for  trapezoid  weighting  function  for  inverse  check 

DO    790    I-1,N1 

FI-I 

IF    (FI.LT. SI)    GO    TO    770 

IF    (FI.LT.S2)    GO    TO    780 

U(  I  )-0.0 

GO    TO    790 

U(  I  )M  .0 

GO    TO    790 

U(  I  )-(  S2-FI  )/(  S2-S1  ) 

CONTINUE 

W(  1  )-U(  1  ) 

GO    TO    620 

END    OF    CALCULATIONS 
WRITE(  I  PR, 81 0) 

F0FM\T(//20X,  10HEND    OF    RUN) 
STOP 
END 


TRANSFORMATION    TO    TRAPEZOID    RULE    APPROX IMATION. CONSTRUCTION    OF    D    MATRIX 
FUNCTION    D(  I,  J  ) 
IF(  J.GE.I*2  )    GO    TO    2 
D-0.5 
RETURN 

D-0.0 
RETURN 

END 


FUNCTION    F(  I,  J  ) 
C0MM0N/BL0CK1 /A, N, I  START, NN, MD, MM 

IF(  I.GE.N-1  )    GO    TO    3 

IF( J.EQ.I )    GO    TO    1 

IF(MOD( J-I,2).EO.O)    GO    TO    2 

IF( J.EQ.N)    GO    TO    5 
F--2./3.*B<  I, J  ) 

RETURN 
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1  F-7./3.*B(  I.J ) 

RETURN 

2  IF(J.EO.N-l)    G0    TO    5 
IF(  J.EQ.N)    GO    TO    4 

F-8./3.*B(  I,J) 
RETURN 

3  F-.5*B(  I  ,J  ) 
KETURN 

4  F-1./3.*B(  I.J  ) 
RETURN 

5  F-2./3.*B(  I, J  ) 
RETURN 

END 
C 
C 
C        COMPUTES  B  MATRIX 

FUNCTION  B  (  I,  J  ) 
C 
C       »*#*«*  **«■*#**#  #*#»*«■*#«■«•*•«•#***  **********  *******  ##«««««#«»»«#««««  »* 

c 

C0MM0N/BL0CK1 /A, N, I  START, NN, MD,  MM 

DATA  PI /3. 1415926/ 
S-I 

T-J 

G0  TO  (10,20,20,30),  MD 
C        B  MATRIX  FOR  CONSTANT  WEIGHTING  FUNCTION 
10     B*PHK2.0*T-A)/2.0 

RETURN 
C        B  MATRIX  FOR  GAUSSIAN  AND  DATA  WEIGHTING  FUNCTION 
20      R-SORT(  (T-A*l  .0  )**2-(  S-l  .0  )*»2  ) 

CALL  WGHT23(  XKI,TKS,R) 

B»(  XKI  *TKS)/(2.0-A) 

R«SQRT(  (T-l  .  0)#*2-(  S-l  .0  )*#2  ) 

IF  (  J.EQ.  I  )  RETURN 

CALL  WGHT23(  XKI,TKS,R) 

B  »B  -(  XKI  *TKS  )/(  2  .  0  -A  ) 

RETURN 
C        B  MATRIX  FOR  TRAPEZOID  WEIGHTING  FUNCTION 
30      R-SORT(  (T-A*l  .0  )**2-(  S-l  .0  )**2  ) 

CALL    WGHT4(  Bl  ,  R  ) 

B"(  2.0*A-1  ,0)*B1 

fi-SORT(  (T-l  .0  )«*2-(  S-l  .0  )**2  ) 

IF    ( J.EQ. I )     RETURN 

CALL    WGHT4(B1,R ) 

B-B-(  2.0*A-1  .0  )*B1 

RETURN 

END 
C 
C  #*#*«****#*###*****#****#»#*###*##*»###***#****#***»##********•»**» 

SUBROUTINE  WGHT23( XKI ,TKS, R ) 
C 
C 

C        CALCULATE  B(I,J)  FOR  DATA  WEIGHTING  FUNCTION 

C       »**«**##**#*##»#*#*##»*#****##**#***»****#**###***#***#«#*###«■*#** 
C 

DIMENSION  TK(400) 

COMMON/BLOCK 1 /A, N, I  START, NN, MD, MM 

C0MM0N/BL0CK2/RN(  4  00  ),W(400  ),V(  400  ),C(  40  0  ),DL,  G 

COMMON    /BL0CK3/    KB,KL 

DATA    PI/3.1415926/ 
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KC»KB*1 

DO  40  K-KC,KL 

IF  (R.LT.  RN(K))  GO  TO  50 
40     CONTINUE 

WF»W(  KL  )♦(  W(  K.L  )-W(  KL-1  )  )*(  R-RN(  KL)  )/(  RN(  KL  )-RN(  KL- 1  )  ) 

LA-KL 

GO  TO  6  0 
50      LA-K-1 

WF-W(LA  )♦(  W(K)-W(LA)  )*(R-RN(LA)  )/(  RN(K)-RN(  LA)  ) 
60  XK0»R#R-RN<  LA )*#2 

IF    (XKO.NE.O.)    GO    TO    70 

XKI-PI#R*R*0.2  5  0*(  W(  LA)*W(LA-1  )-(  RN(  LA  )*RN(  LA- 1  >  )*(  W(  LA  )-W(  LA-1  )  )/ 
2(  (  RN(  LA  )-RN(  LA-1  )  )  )  ) 

GO    TO    8  0 
7  0  XK1  «SQRT(  XKO  ) 

XK2"RN(  LA)/XK1 

XK3  -(  R*R*ATAN(  XK2  )*RN(  LA  )*XK1  )*RN(  LA  )*XK1*»3/1  .  5 

XK4«<  W(  LA-1  )-(  DL/(  R-RN(LA)  )  )*(  W(  LA  )-WF)-W(  LA)  )#XK3/DL 

IF     (LA.EQ.l)     XK4-0.0 

XK5«PI*R*R*0.250#(  WF*W(  LA)-(  R*  RN(  LA  )  )/(  R-RN(  LA  )  )#(  WF-W(  LA  )  )  ) 

XKI «XK4*XK5 
80  TKS-(  W(  2)-W<  1  )  )*R*R*R/(  1  ,  5»DL  ) 

IF    (LA.EO.l)    TKS-C  WF -W(  1  )  )*R»R*R/(  1 .  5*DL  ) 

KK-LA-1 

IF    (KK.LT.KC)    RETURN 

DO    110    K«KC,KK 

TKD-R«R-RN<  K  )**2 

IF    (TKD.NE.O,)    G0    T6    90 

TKA-PI*R*R*0.5*RN(  K  )/DL 

GO    T0    1 00 
90  TKA-RN(  K)*R*R*ATAN(  RN(  K  )/SQRT(  TKD  )  )*RN(  K)*RN(K  )*SQRT(  TKD  )♦<  SQRT(  TK 

2D**3) )/l .5 
100  TK(  K)«(  \W  K*l  )*W<  K-l  )-2.0*W(K)  )#TKA/DL 

TKS-TKS*TK(K  ) 
110  CONTINUE 

RETURN 

END 
C 
c  «##**»»***##****##***»#*#*#»*****»«*#*#«*****#*####*#****#****#*** 

c 

SUBROUTINE    WGHT4(B1,R) 
C 

COMMON    /BL0CK5/S1.S2 

DATA    PI/3.1415926/ 
C  CALCULATES    B(  I , J )    FOR     TRAPEZOID    WEIGHTING    FUNCTION 

C 

IF    (R.LE.S1)    GO    TO     130 

IF    (R.LE.S2)    GO    TO    120 

B2-S2*SQRT( R»R-S2*S2 ) 

B3"S1*SQRT(  R»R-S1*S1  ) 

B6-S1*S1/B3 

B7-S2*S2/B2 

B4»R#R»ATAN(  B7  ) 

B5-R*R*ATAN(  B6  ) 

IF(  SI .LT. 0.000001  )    GO    TO    115 

Bl-(  B3*B5)*0.5*(  B2-B3»B4-B5  )*0.  5*S2/(  S2-S1  )♦((  B2/S2  )**3-(  B3/S1  )##3 
2)/(  3.0*(S2-S1  )  ) 

RETURN 
115  B1-(B2*B4  )*0.5*((  B2/S2  )**3-R**3  )/(  3.0*S2  ) 
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RETURN 
120  IF(  SI  .LT. 0.000001  )    G6    T0    140 

B3-S1#SQRT(  R#R-S1*S1  ) 

B6-S1*S1/B3 

B5-R*R*ATAN(  B6  ) 

Bl*(  B3*B5  )*0.5*(  -B3-B5*R#R* PI/2.0  )*S2*0.5/(  S2-S1  )-(  (  B3/S1  )#*3  )/(  3  . 
20*(  S2-S1  )  ) 

RETURN 
130  Bl "R*R*PI*0.25 

RETURN 
140  Bl  "R*R*PI*0.25-(  R#*3  )/(  3.0*S2  ) 

RETURN 

END 
C 
C 

SUBR6UTINE    CGEN    (I) 
C  **«*****#****#****#**#*«•*##*#*#*#**#*#*»*#***#*»#**##**«#*##*****» 

C  CALCULATE    ELEMENTS    6F     ITH    ROW    0F    MATRIX    C 

C 

C0MM0N/BL«CK1/A,N,ISTART,NN,MD,MM 

C0MMON/BLOCK2/RN(  400  ),W(  400  ),V(  400  ),C(  400  ),DL,G 

DATA    PI/3.1415926/ 

S-I 

SN-N 
C  CHECK    F0R    LOW    ORDER    APPROXIMATION 

IF    (MM.NE.O)    GQ     TO     140 

IF    (I.EO.N-1)    GO    TO    120 

IF    (I.EQ.N)    GO    TO    130 

DO    110    J-I,N 

U-J 

IF    ( I.EQ. 1 )    GO    TO    50 

IF    (  I-J  )     10,30,40 
10  IF    (  M0D(  J-I,2).EQ.  0  )    GO    TO    20 

IF    (  J.EQ.N)    GO    TO    80 

Ul-J 

U2-J-2 

U3-J-1 

Fl  -SQRT(  Ul  *U1  -(  S-l  .0  )*#2  ) 

F2»S0RT(  U2*U2-<  S-l  .  0  )**2  ) 

T-(  Fl**3-F2**3  )/3.0*( S-l .0  )**2*( Fl -F2)-U3*( U1*F1-U2*F2 )-U3*(  S-l . 0) 
2**2*AL0G(  (  Ul  *F1  )/(  U2+F2  )  )*U3*U3*(  F1-F2  ) 

V(  J)*F1  -F2-T 

GO    TO    110 
20  Ul-J-1 

U2-J-3 

U3-J-2 

Fl  *SQRT(  Ul  #U1  -(  S-l  .  0  )**2  ) 

F2-SQRT(  U2*U2-(  S-l  .0  )**2  ) 

T"(  Fl**3-F2*#3  )/3. 0*( S-1.0  )**2*( Fl -F2  ) -U3#( U1*F1- U2*F2 ) -U3*( S-l  . 0) 
2**2*AL0G(  (  Ul  *F1  )/(  U2*F2  )  )*U3#U3*(  Fl  -F2  ) 

D-0.5#(  U1*F1  -U2*F2M  S-1.0  )**2*AL0G(  (  U1*F1  )/(  U2*F2  )  )  )-U3*(  Fl -F2  ) 

V(  J)"0.5*(  D*T) 

IF     (J.EQ.N)    GO    TO    110 

IF    (J.EQ.N-1)    GO    TO    90 

Ul-J*l 

U2-J-1 

U3-J 

Fl  -SQRT(  Ul  *U1  -(  S-l  .0  )**2  ) 

F2"S0RT(  U2*U2-(  S-l  .  0  )**2  ) 

T-(  Fl»*3-F2**3  )/3.0*( S-l .0  )**2*  (  Fl -F2  ) -U3*<  U  1  »F1-U2*F2  ) -U3#(  S-1.0) 
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2#*2«AL0G(  (U1*F1  )/(  U2*F2)  )*U3»U3*(  F1-F2) 

D-0.5*(  U1*F1  -U2*F2*(  S-l  .  0)**2*AL0G(  (  Ul  *F1  )/(  U2  *F2  >  )  )-U3  *<  Fl  -F2  ) 

V(  J  )-V(  J  )-0.5*<D-T  ) 

GO    TO     1 10 
30  Ul-I*l 

U2-I-1 

U3"I 

Fl  »SQRT(4.0*S) 

T«Fl##3/3.0*(  S-l .0  )**2*F1-U3*U1*F1-U3*( S-l. 0  )**2*AL0G(  (  U1*F1 )/U2 ) 
1  ♦U3#U3*F1 

D-0.5*(  U1#F1  ♦(  S-1.0  )#*2*AL0G(  (  U1*F1  )/U2)  )-U3*Fl 

V(  J  )--0.5*(  D-T) 

GO    T0    110 
40  V(J)»0.0 

G0    TO    110 
50  IF    (I-J)    60,70,100 

60  V<  J  )  - 1  .  0 

GO    TO    110 
70  V(  1  )»0.5 

GO    TO    110 
8  0  F1»SQRT(<  D-l  .0  )*»2 -(  S-1.0  )**2  ) 

F2-SQRTC  (  U-2  .0  )*#2  -(  S-1.0  )**2  ) 

V(  J)-(  F1-F2  )*0.5 

G0    TO    110 
90  Fl  -SQRT(  U*U-<  S-1.0  )**2  ) 

F2-SQRT(  (  U-l  .0  )»*2-(  S-l  .  0  )##2  ) 

V(  J  )  -V(  J  )  ♦(  Fl  -F2  )*0  .  5 

GO    TO    110 
100  V(J)*0.0 

110  CONTINUE 

RETURN 
C  END    POINTS    CALCULATIONS 

120  F4-SQRT(2.0*(  SN-1.0)-1  .0  ) 

V(  N-l  )-0.5*F4 

V(  N)-0.5*F4 

RETURN 
C        LOW-ORDER  APPROXIMATION 
130    V(  N)-0.5*S0RT(  2.0*SN-1.0  ) 

RETURN 
140    DO  190  J-I  ,N 

U-J 

IF     (  I-J  )    150,160,170 
150  IF    (J.EQ.N)    GO    TO    180 

V(  J)»(  SQRT(  U*U-(  S-1.0  )*#2  )-SQRT(  (  U-  2.  0  )**2-(  S-  1  .  0  )**2  )  )  *0.  5 

GO  TO  190 
160    V(  J)»(  SQRT(  2.0»S-1  .0  )  )*0.5 

GO  TO  190 
170    V(J)»0.0 

GO  TO  190 
180    V(  J)-0.5*(  SQRT(  (  U-l  .0  )#»2-(  S-1.0  )**2  )-SQRT(  (  U-2.0  )#*2-(  S- 1  .  0  )**2  )  ) 
190    CONTINUE 

RETURN 

END 
C 

C     CORRECTION  FOR  THE  EFFECTS  OF  THE  WIDTH  USING  THE  METHOD  OF  TAYLOR 
C     AND  SCHMIDT 

SUBROUTINE  WIDTH(  FF,  LW,  WQ,  KO  ) 

DIMENSION  FF(4  00  ),F(  400  ),LW(  10  ),  WQ<  10  ) 

COMM0N/BLOCK1 /A, N, I  START, NN,  MD,  MM 

AA-756. 
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LW(K0*1  )-N*l 

DG    1     K-1.K6 

NW"LW(  K  )*2 

KW-LW(  K*l  )-4 

D6    305    J-    NW.KW 
305    F(J)     -    FF(J)-WQ*    (FF(  J-2  )#1  20  . -FF<  J- 1  )»87.-FF<  J)»96.-FF<  J*l  )*12.* 
1     FF(  J*2  )*60.*FF(  J*3  )*15.   )    /    AA 

J0-LW(K  > 

IMAX-LW(K*1  )-l 

F(  JO)     -    FF(  J0)-WQ»(3  30.*FF(  JO  ) -381  .  #FF(  JO*l  )-26  4.#FF(  JO  *2.  )♦  1  56.  *FF 
1( J0*3  )*354.»FF(  J0*4  )-195.*FF< J0*5)  )/AA 

F(  J0*1  )    -    FF(  JO*l  )-WQ»(225.*FF(  JO  )-234.*FF(  J0*1  )-180.*FF(  J0*2)*72. 
1*FF(  J0*3)*207.*FF(  J0*4  )-90  ,*FF(  JO*  5  )  )/AA 

F(IMAX-2)     -    FF    (IMAX-2)    -    WQ»(     1  5.  *FF(  IMAX- 5  )♦     60  .  *FF(  I  MAX-  4  )  -1  2  .  # 
lFF(IMAX-3)-    96.  *FF(  IMAX-2  )-    87 . »FF(  IMAX- 1  )♦ 1 20. *FF(  IMAX  )  )    /AA 

F(IMAX-l)     -    FF    (  IMAX-1  )     -    WQ#( -90.  »FF(  IMAX-5  )♦  207.  *FF(  I  MAX-4  )  *72  .  * 
1FF(  IMAX-3 )-l 80.*FF(  IMAX-2 )-234.*FF(  IMAX- 1  )*225.*FF(  IMAX  ))    /AA 

F(IMAX)        ■    FF(IMAX)       -WQ*( -1 95 . *FF(  IMAX-5 )    ♦ 354. *FF(  IMA X-4 )     *156.» 
1FF(  IMAX-3)     -264.  *FF(  IMAX-2  )    -381  .  *FF(  IMAX-  1  )*330  .*FF(  IMAX  )  )    /AA 

CONTINUE 

D0    2    J*1,N 

FF(  J)-F(  J) 

RETURN 

END 


CHAPTER  6 

1.   TYPICAL  CALCULATIONS  AND  DISCUSSION  OF  VARIOUS  ALGORITHMS 

In  Figures  1-4,  unsmeared  data  calculated  from  the  measurements 
supplied  by  J.  J.  Weeks  of  this  laboratory  are  presented.   The  measure- 
ments were  performed  on  melt-crystallized  polyethylene. 

In  Figure  1,  results  based  on  very  low  scattering  angles  are  shown. 
Results  based  on  Nystrom  algorithm  (forward  interpolation)  are  shown  by 
a  solid  line.   Results  based  on  the  rectangle  rule  approximation  and  on 
the  trapezoid  rule  approximation  (uncorrected  for  perturbation  effects) 
are  shown  by  sets  of  points.   As  is  evident,'  even-odd  oscillations  start 
to  be  noticeable  at  larger  angles. 

In  Figure  2,  results  for  larger  angles  are  shown  (20>.OO8  rad.). 
The  bottom  set  of  points  represent  the  input  data.   The  desmeared  data 
are  shown  for  two  algorithms:   rectangle  rule  approximation  and  Nystrom 
algorithm.   While  both  methods  of  interpolations  fall  on  practically  the 
same  curve,  the  Nystrom  algorithm  seems  to  lead  to  more  numerically 
stable  results  (smaller  oscillations  from  the  smooth  curve) . 

As  we  pointed  above,  the  results  for  desmearing  data  are  very 

sensitive  to  the  numerical  perturbations  in  the  input  data.   If  the  data 

show  excessive  fluctuations  (resulting  from  experimental  error  or  random 

noise  effects),  those  will  be  magnified  when  desmeared.   To  demonstrate 

this,  we  present,  in  Figure  3,  results  obtained  from  experimental  data 

with  certain  amounts  of  random  noise  (bottom  curve  in  Figure  3) .   The 

upper  curve  and  points  show  results  of  desmearing  data  represented  by 

the  bottom  curve.   The  sharp  fluctuations  in  the  desmeared  data  are 

present,  no  matter  what  algorithm  is  used  in  the  calculations.   In  this 
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figure,  the  solid  line  represents  the  results  based  on  the  Nystrom 
algorithm.   These  results  are  somehow  more  smooth  than  the  ones  based  on 
other  algorithms.   On  the  other  hand,  the  trapezoid  rule  approximation, 
even  after  correction  based  on  perturbation  method  (smoothing  the  results) 
were  applied,  lead  to  the  biggest  scatter  in  the  calculated  data. 

In  Figures  1-3  a  constant  weighting  function  was  assumed.   In  Figure 
4,  trapezoid  weighting  function  was  assumed.   As  expected,  the  presence 
of  a  weighting  function  has  an  effect  on  stabilizing  the  numerical 
solution.   The  unsmeared  data  are  devoid  of  scatter,  and  all  three 
interpolating  methods  employed  (see  captions  to  Figure  5  for  their 
description)  fall  into  identical  curve. 

In  Figure  5  we  demonstrate  the  effect  of  desmearing  on  peak  formation 
in  the  scattering  curve.   The  original  (input)  data  were  provided  by  D. 
Mclntyre  of  the  Polymer  Institute,  University  of  Akron.   The  measurements 
were  performed  on  a  triblock  polymer  of  polystyrene-polybutadiene- 
polystyrene. 

The  lower  curve  in  Figure  5  represents  the  original  data.   The 
upper  curve  shows  the  desmeared  data,  obtained  with  the  Nystrom  algorithm 
(forward  interpolation) .   Of  interest  is  the  existence  of  the  first  peak 
in  the  desmeared  curve  (indicated  by  an  arrow) ;  there  is  no  corresponding 
maximum  in  the  input  data,  only  an  inflection  point.   We  also  notice  the 
following  facts:   1)  the  maximum  in  the  input  data  are  greatly  magnified 
when  these  data  are  desmeared;  and  2)  the  desmearing  of  the  scattering 
curve  tends  to  shift  the  locations  of  the  peaks  toward  slightly  larger 
angles. 
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In  all  these  Figures,  only  a  part  of  the  scattering  curves  are 
shown.   In  Table  3  and  4  typical  calculations  are  presented,  together 
with  the  input  data.   Fragments  of  the  input  and  output  data  presented 
in  Table  3  and  4  are  depicted  in  Figures  4  and  5. 

2.   TEST  FUNCTIONS 

Two  functions  were  employed  using  a  constant  weighting  function: 
1)  I(x)  is  Gaussian.  This  function  was  employed  in  order  to  compare 
the  accuracy  of  the  various  algorithms  employed  and  in  order  to  check 
the  effect  of  the  cut-off  on  the  relative  error.   The  function  employed 


is: 


y  -i2/1250  .    .  n   ,        im 

I.  =  e         for  i=0,  1,  ....  100, 

i 


(49) 


Tabulated  values  of  I. /I.,  using  various  algorithms,  are  shown  in  Table 
5.   The  exact  results  are  computed  from: 


—  =  0.01595 

I. 

i 


L    exp[-(10,000-i2)/1250l) 

;  /io,ooo-in3/2  ^2l 

\     1250   J 


(50) 


As  is  evident,  Nystrom  algorithm  leads  to  best  agreement  with  exact 

results.   At  x=90 ,  the  error  because  of  the  cut-off  at  x=100,  is  still 

about  1%. 

2)   The  function  was  taken  from  P.  Schmidt  [10].   It  represents  the 

infinite-slit  smeared  scattering  function  for  monodisperse  spheres. 

This  function  corresponds  to  the  perfect  collimation  intensity  function 

with  a  rapid  succession  of  maxima  of  diminishing  intensities.   Hence,  it 

provides  a  good  test  for  the  numerical  method  employed. 

The  perfect  collimation  intensity  I(x)  for  spheres  is  given  by 

I.  =  9  [sin(ia)  -  ia  cos(ia)]2/,.   6  (51) 

i  via,} 
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i  is  a  scattering  angle.   a=0.25;  i=3,4,  ...,  80.   The  results  are  shown 
in  Figure  6.   In  this  figure,  ia=l  corresponds  to  a  scattering  angle  of 
4  miliradians.   The  solid  line  represents  eq.  (51).   The  points  are  the 
calculated  intensities.   Method  of  forward  interpolation  (Nystrom 
algorithm)  was  employed  in  the  calculations. 

As  is  evident  from  Figure  6,  there  is  a  close  agreement  between  the 
calculated  and  exact  results,  almost  up  to  i=76.   The  cut-off  of  calcula- 
tions was  at  i=80  (corresponding  to  scattering  angle  of  80  miliradians) . 


The  author  wishes  to  thank  Dr.  John  Smith  for  initiating  the  computer 
program  and  to  Mr.  James  Weeks  for  helpful  suggestions  about  handling 
experimental  information  and  also  for  providing  experimental  data.   Mr. 
Martin  Cordes '  help  in  preparation  of  this  Technical  Note,  and  in  particu- 
lar, his  help  with  the  analysis  of  the  mathematical  background  to  this 
work,  is  very  appreciated. 
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TABLE  1 


SYMBOL 


FORMAT 


CONDITION 


COMMENT 


NN 

ISTART 
IENTR 
NSUB 


4110 


First  card. 


IENTRO  indicates 
that  there  is  no 
interpolation  of 
input  data. 


LW* 

WM2> 


110,  F10.5 
(per  card) 


Number  of  cards  - 
number  of  slit- 
width  corrections 
+  1  (blank) 


Last  card  must  be 
blank.   If  no  slit- 
width  correction 
applies  a  single 
blank  card  must  be 
inserted. 


MD 

IDC 

INV 

MM 


4110 


Single  card 


G 

F10.2 

MD  =  2** 

RN* 

w* 

2F10.2 

MD  =  3** 

Series  of  cards. 
Last  card  must  be 
blank. 

SI 
S2 

2F10.2 

MD  =  4** 

XI* 

8F10.0 
card) 

(per 

Last  card  must  be 
blank. 

Dimensioned  variables. 

These  cards  are  inserted  only  according  to  the  directive  given  by  the 
value  of  the  integer  MD . 
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Table  2 

ARRAYS  USED  IN  PROGRAM 

SYMBOL  DESCRIPTION 

XI  Input  Vector 

YY  Output  Vector,  rectangular  rule  approximation 

YI  "      "    ,  trapezoid  rule  approximation 

YK  "      "    ,  smoothing  correction  applied 

to  YI 

YL  ,  smoothing  correction  applied 

to  YK 

YM  "      "    j  Nystrom  algorithm  (forward 

interpolation) 

XX  Temporary  input  vector,  for  data  interpola- 

tion. 

ZI  Input  vector,  recalculated  by  solution 

checking  (used  only  if  INV^O) 

W  Input  weighting  function 

RN  Input  coordinates  of  W 

U  Temporary  storage  for  weighting  function 

for  solution  checking  (used  only  if 
INV^O) 

V  Array  of  y . .    for  given  i 

N        XJ 

C  /j      V. .  X. 

3=1 
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SPECIFICATIONS   • 
NrTOT/'L  RANGE  OF  DATAr=10C     NN» NUMBER  OF  PAT*  RFAD=!00     ISTApTr STARTING  POTNT»=   0 

If.  IS  TRAPEZOID'Sl=    25.000#S2=     60. GOO 


C-VALUE.S.. 


0 

.984  02+02 

42 

.29645+^1  . 

64 

.92378+00 

1 

.83610+0? 

43 

.28816+01 

85 

.88917+00 

2 

.76661+02 

44 

.27,:»4  9  +  nl 

66 

.35343+00 

3 

.70560+02 

45 

.27083+01 

87 

.81716+00 

4 

.64307+02 

46 

.26247+01 

68 

.77047+00 

5 

.5'7646+02 

47 

.254  63+01 

69 

.7414^+00 

6 

.50P16+0? 

48 

.24733+01 

9C 

.70166+00 

7 

.43717+02 

49 

.24040+01 

91 

.66001+00 

8 

.37058+02 

50 

.23397+01 

92 

.61660+00 

.  9 

.30766+02 

51 

.22750+01 

93 

.5701 3+0 0 

10 

.25355+0? 

5? 

.22141+OJ 

94 

.5197o+00 

11 

.2101.4  +  02 

53 

.21534+01 

95 

.46425+00 

12 

.17652+02 

54 

.20961+01 

9b 

.40143+00 

1? 

.15193+0? 

55 

,2044'++01 

97 

.32734+00 

1'+ 

.13403+02 

56 

.1994!;  + 0.1 

98 

.23089+00 

15 

.12239+02 

57 

.19471+01 

99 

. 1156^  +  00 

16 

. 1] 369+02 

56 

.19016+01 

17 

.10643+02 

59 

.18567+01, 

18 

.99876+01 

60 

.18166+01 

19 

.94170+01 

61 

.1777'.  +  "! 

20 

, 864  60+ :H 

62 

.17394+01 

21 

.8342-3  +  01 

63 

.17010+01 

22 

.772^5+01 

64 

.16647+01 

23   . 

.72106+01 

65 

.16265+01 

24 

.67074+01 

66 

.15fi6"7  +  0! 

25 

.62452+0.1 

67 

. 15471+01 

26 

.5i!335  +  01 

68 

.150P6+01 

27 

.54593+01 

b° 

.14685+01 

2f 

.51182+01 

70 

. 14290+01 

29 

.48080+ri.l 

71 

.  13897+01 

30 

.45425+01 

72 

.13514+01 

31 

.43229+01 

73 

.13136+01 

32 

.41454+01 

74 

.12755+01 

33 

. 39885+0 1 

75 

.12373+01 

34 

.38337+03 

7o 

.11995+01 

35 

.3b922+o'l 

77 

.11635+01 

36 

.35697+01 

78 

.112.84+01 

37 

.34^19+01 

79 

.in'334+ni 

36 

.33454+01 

BO 

.10592+01 

39 

.32309+01 

81 

.1025^+01 

40 

.31417+01 

82 

.99ORS+00 

41 

.30501+01 

63  ~ 

.95867+0  0 

62 


ARBITRARY  UNITS 

INPUT  DAT/*  EXTRAPOLATED  TO  ZERO  ANGLE 


INPUT  VFCTOR   2F»OTH  ORDER   TRAPEZOID    TRAPE70ID+PERT.  "T? APF/.0ID+2PFPT.  NYSTROM  INTER. 


y  i  ( i ) 


YY(I) 


Yim 


YK(I) 


YL(T) 


YM  ( I ) 


0 

374000+02 

.Ql+l^crn-Ol 

.166500+02 

.168500+02 

.168500+02 

.000000 

1 

138P00+C2 

.14-7069+01 

.1986  02+01 

.42.916^+01 

.584756+01 

,326784+01 

C             * 

1C3P00+02 

.779972+00 

.9S5_756+00 

.106887+01 

.153638+01 

.102463+01 

•7 

868000+01 

.569320+00 

.604! 59+00 

.661160+00 

.69978,2+00 

.656559+00 

4 

785C0O+01 

.4-71755+00 

.534451+00 

.525176+00 

.529571+00 

.515821+00 

5 

667000+01 

.395878+00 

.409059+00 

.425564+00 

.433837+00 

. 4*3296+00 

6      « 

581 0  0  0+01 

.348198+0.0 

.382697+00 

.375591+00 

.373827+00 

.372754+00 

7 

4  890  0  0+01 

.283215+00 

.313699+tJO 

.3.15037+00 

.315518+00 

.315484+00 

fc      < 

41.10000+01 

.236164+00 

.252731+00 

.257370+00 

.25^510+00 

.260375+00 

9 

320000+01 

.181866.+00 

.21 9598+00 

.? 12543+00 

.21.0040+00 

.207721+00 

10 

.240000+01 

.132104+00 

.144134+00 

.162701+00 

.156433+00 

.155251+00 

13 

1*1000+01 

.93578"- 01 

. 1 20075+00 

.115253+00 

.114719+00 

.111273+00 

12 

1320  011  +  01 

,631338-Q] 

.670  825-01 

,74598,4-01 

.781095-01 

.770785-01 

13 

990000+00 

.(+27127-01 

.591851-01 

.580106-01 

.543833-01 

.516637-01 

11 

74-0000+00 

,260245-Cl 

.262403-01 

.316592-01 

.341725-01 

.331945-01 

15 

595000+00 

.163302-01 

•258007-01 

,233878-01 

.232936-01 

.216987-01 

If 

520000+00 

.144720-01 

.  1U8516-01 

.146513-01 

.160945-01 

.162414-01 

17 

4-65000  +  00 

.12357°-P1 

.18092.4-01. 

.14-9741-01 

.140311-01 

.132851-01 

ie 

420000+00 

.102354-01 

.66231.7-02 

.973880-02 

.10*674-01 

.111945-01 

is. 

390000+00 

.974 181-02 

.138477-01 

.112750-01 

.104518-01 

.985671-02 

2(J 

359000+00 

.823824-02 

.563590-02 

,7«71c'P-0? 

.871470-02 

.904751-02 

21 

341000+00 

.949694-02 

.1O34C6-01 

.952526-02 

.9071.50-02 

,89988°-02 

?2 

314000+00 

.771770-02 

.815329-02 

.0455^8-02 

.845811-02 

.054245-02 

£3 

273000+00 

.71Q86°-o? 

.728211-02 

.739950-02 

.749181-0? 

.748865-02 

24 

255000+00 

.637677-02 

.711527-02 

.689691-02 

.680443-02 

.677101-02 

2b 

2.20C00+00 

.550336-0? 

. 563 p 26-0 2 

.583946-02 

.591585-02 

.5Ql?50-02 

26 

201000+00 

.464911-02 

.536845-02 

.524031-02 

.521537-02 

.516678-02 

27      < 

177000+00 

.429566-02 

,a?2Q76-02 

.449151-02 

.455437-02 

.458257-02 

28      < 

160000+00 

.360249-02 

.426156-0? 

.41 19^0-0? 

.406970-02 

,405°75-02 

29 

140COO+00 

.319542-02 

.334341-05" 

.344711-0? 

.347680-02 

.3^8*36-02 

30 

1240  00+00 

.261458-02 

,304744-02 

.pc.5248-02 

.292646-02 

.287215-02 

31 

11 0000+00 

,210894-02 

.215172-02 

.230175-02 

.235937-02 

.232659-02 

32 

leipoo+oo 

.164520-32 

..''0361  7-0? 

.199677-02 

.200765-02 

.107443-02 

33 

9500  00-01 

. 17716^-02 

.165423-02 

.175709-02 

.179647-02 

.181919-02 

34 

690000-01 

.159073-02 

.  18.894-3-02 

.175366-0? 

.170054-02 

.167821-02 

35 

520000-01 

.138472-02 

.131002-02 

.143149-02 

.147628-02 

.148201-02 

36 

775000-01 

.129005-02 

. l45°4]-0? 

.137606-02 

.135751-02 

.1^3433-02 

37 

735000-01 

.119334-02 

.112068-02 

. 120136-02 

.123090-02 

.124*21-02 

3c       < 

700000-01 

. 114633-02 

.126600-02 

.1201*9-02 

.118016-02 

.117199-02 

39      < 

6-65000-01 

.104-308-02 

■ .102666-02 

.107202-02 

.108583-02 

.109211-02 

4  0 

630000-01 

.972527-0  3 

.1 05949-02 

.1 02503-02 

.101549-02 

.100420-02 

"1 

6050  00-01 

.904241-03 

.865560-03 

."?0776-03 

.934114-03 

.931*36-03 

42 

580C  00-01 

.862476-03 

.922921-03 

.896546-03 

.892015-03 

.88272?-03 

43       1 

565( 0  0-01 

.670865-03 

.802034-03 

."45126-03 

.860664-03 

.871119-03 

44 

545000-01 

.8502.42-03 

.939697-03 

.886935-03 

.866692-03 

.86370^-03 

4b 

520000-01 

.810829-03 

.760788-03 

.8ri7286-03 

.820963-03 

.831494-03 

4b 

495C0C-01 

.757400-03 

.860871-.03 

.809700-03 

.791928-03 

.783*44-03 

47 

470000-01 

.70510+-03 

.653930-03 

.705478-03 

.723491-03 

.730357-03 

48 

4-50000-01 

.666619-03 

.756278-03 

.709334-03 

.694042-03 

.6*3863-03 

49 

.430000-01 

.620734-03 

.576961-03 

.6?1438-03 

.637783-03 

.642580-03 

63 


JOT- 


50 

.^15000-01 

.609593-03 

.664^07-03 

.631612-03 

.621426-03 

.615614-03 

51 

.400000-01 

.577740-03 

.5*4679-03 

.580671-03 

.589174-03 

.594325-03 

52 

.384000-01 

.56370  4-0^ 

.600800-0"* 

.^B07H3-03 

.574135-03 

.572.98q-03 

53 

.370000-01 

.533132-03 

.526609-0"! 

. ^41148-03 

.544431-03 

.547322-03 

54 

.352000-01 

,49095"*-03 

.5.39654-0-* 

,<^Pl2'J6-03 

.515779-03 

,510«65-03 

55 

.340000-01 

.472415-03 

.4'-l-2?52-03 

."68540-03 

.478327-03 

.4*1280-03 

56 

.329000-01 

.44866:>-03 

.^02579-03 

,4745c:4>-03 

.465284-03 

.460255-03 

57 

.-318000-01 

.431413-03 

.394756-03 

.4249"6-03 

.435754-03 

.441402-03 

56 

.310000-01 

.4210b9-03 

.468067-0  3 

.440183-03 

.429970-03 

.426153-03 

59 

.299000-01 

.387283-03 

.374071-03 

.394141-03 

.400785-03 

.402^00-03 

60 

.291000-01 

.37611,0-03 

.400496-03 

.3«7963-03 

.384961-03 

.381785-03 

61 

.286000-01 

.363479-03 

.351723-03 

."*63770-03 

.367908-03 

.369300-03 

62 

.280000-01 

.354498-03 

.375234-03 

.3644  04-03 

.361029-03 

.358317-03 

63 

.276000-01 

.346108-03 

.333761-03 

."44789-03 

.348922-03 

.349680-03 

64 

.272000-01 

.345351-03 

.358455-03 

.349971-03 

.347699-03 

.346750-03 

65 

.269000-01 

.346963-03 

.332246-03 

.341520-03 

.344248-03 

.346888-03 

6b 

.26'i000-01 

.339672-03 

.361680-03 

.349438-03 

.344647-03 

,34231q-03 

67 

.258000-01 

.323494-03 

.317664-03 

."23610-03 

.332783-03 

.334301-03 

6a 

.254000-01 

.330633-03 

.339323-03 

.332817-03 

.330728-03 

.330437-03 

69 

.249000-01 

.320897-03 

.321944-03 

.3244°2-03 

.324945-03 

.325892-03 

70 

.243000-01 

.315007-03 

•319P51-03 

."18885-03 

.318525-03 

.318766-03 

71 

.238000-01 

.305^84-03 

.311964-0^ 

.311113-03 

.310696-03 

.310042-03 

72 

.232000-01 

.297^71-1^ 

.249005-03 

.300820-03 

.301803-03 

.301022-03 

73 

•228000-01 

.294354-03 

.296^37-03 

. 2964^1-03 

.296437-03 

.296411-03 

74 

.223000-01 

.290231-03 

.291771-03 

.202110-0"* 

.292062-03 

.293161-03 

7b 

.210,000-01 

.28349^-0^ 

•2886Ql-03 

.287473-03 

.286657-03 

.286917-03 

76 

.212000-01 

.271°85-0  3 

.278305-03 

,'77q?q-03 

.277754-03 

.278203-03 

77 

.207000-01 

.264347-03 

.265664-03 

.'6,73^2-03 

.268239-03 

.268631-03 

76 

.20^000-01 

.259156-03 

.263030-03 

.?6217-'-03 

.261903-03 

.264378-03 

79 

.1990  0  0-01 

.251700-03 

.256282-03 

.2553"'9-03 

.255137-03 

.254322-03 

80 

,195000-01 

.2415bl-03 

.248] 18-03 

.247126-03 

.246824-03 

.245810-03 

81 

. 192000-01 

.234621-03 

.235003-03 

.237062-03 

.238189-03 

.237057-03 

82 

.'190000-01 

,23242'-i-0  3 

.234239-03 

.233763-03 

.2337^,0-03 

.233632-03 

83 

.188000-01 

.228652-03 

.2"06]9-O" 

.''30566-03 

.230254-03 

.230467-03 

84 

.185000-01 

.22114U-P3 

.226685-03 

,2?54ou-o3 

.225178-03 

.224443-03 

6E 

,183000-01 

.218829-03 

.215604-03 

.218526-03 

.219658-03 

.220067-03 

86 

.181000-01 

.214170-03 

.222054-0" 

.216351-03 

.217154-03 

.216^*67-03 

87 

.179000-01 

.212  531-03 

.206286-03 

.2109«5-03 

.212501-03 

.23  3285-03 

bt 

.177000-01 

.20670^-03 

.2.18777-03 

.212672-03 

.210677-03 

.209057-03 

89 

.175000-01 

.20577'4-C3 

.194641-03 

.502375-03 

.205210-03 

.20628^-03 

90 

.374000-01 

.204590-03 

'  .2l6°07-03 

.2090Q0-03 

.206415-03 

.204789-03 

91 

.172000-01 

.202395-03 

.102273-03 

.109753-03 

.202387-03 

.203O&0-03 

92 

.1710  0  0-01 

.203755-03 

.2.12517-0" 

.206222.-03 

.204387-03 

.204361-03 

93 

.170000-01 

,2062fc:f:>-03 

.194903-03 

.201679-03 

.204083-03 

.205349-03 

94 

.169000-01 

.210806-03 

.2175°6-03 

.211560-03 

.209984-03 

.200818-03 

95 

.168000-01 

.221114-0  3 

.204026-03 

.21 19^2-03 

.215804-03 

.192538-03 

96 

.167000-01 

.246700-03 

.238202-0"* 

.2353^8-03 

.237600-03 

.186856-03 

97 

.166000-01 

.314886-03 

.255199-03 

.272262-03 

.279718-03 

.343^93-03 

96 

.165000-01 

.372314-03 

.374^73-03 

.353924-03 

.281326-03 

.723^43-03 

99 

.164000-01 

.370056-03 

.370056-03 

.000000 

.000000 

.111279-02 

64 


BACK  CHFCK^ 
INPUT      Rf 


COMPUTED 


XT(I) 


ZI(T) 


0 

.374-00+0-? 

.24612+0? 

50 

.41500-01 

40985-01 

1 

.13800+02 

.1309^+02 

51 

.40000-01 

391.21-01 

2 

.10300+02 

.10145+02 

52 

.3840  0-0  1 

37905-01 

3 

.86800+01 

•  85G4<~  +  01 

53 

.37000-01 

36136-01 

4 

.76500+°1 

.7F566+01 

54 

.352.00-01 

34649-01 

5 

.66700+01 

.66062+01 

55 

.34000-01 

33254-01 

6 

.58100+01 

.58016+01 

56 

.32900-01 

323  3^-01 

7 

.48900+01 

.40465+01 

57 

.31800-01 

31097-01 

8 

.40000+01 

.39053+01 

56 

.31000-01 

30489-01 

o 

.32000+01 

.31603+01 

59 

.29900-01 

29151-01 

10 

.24000+01 

.23713+01 

60 

.2910  0-01 

28603-01 

11 

.18100+01 

.18002+01 

61 

.28600-01 

27Q22-01 

12 

.13200+01 

.13040+01 

62 

.28000-01 

27502-01 

13 

.99000+00 

.9R528+00 

63 

.27600-01 

26951-01 

14 

.74000+00 

.7188a+CC 

64 

.27200-M 

26724-01 

IF. 

.59500+00 

.59830+00 

65 

.26900-0  1 

26321-01 

16 

.520 0  0  + C' 

.50683+00 

66 

.26400-0  1 

25902-01 

37 

.46500+00 

.46964+00 

67 

.25800-01 

25224-01 

la 

.42000+00 

„4<!<jc;,u  +  00 

68 

.25400-jl 

,24954-01 

19 

.39000+0° 

.3953°+00 

69 

.24900-0.1 

.24336-01 

20 

.3590  0+0° 

.34873+00 

70 

.24300-^1 

.23853-01 

21 

.34100+00 

,3a;\4P  +  00 

71 

.23600-01 

.23?6t;-01 

22 

.  U4CO+00 

.30369+00 

7? 

.23200-01 

•2274U-01 

23 

.2730  0+01 

.27738+00 

73 

.22HO°-01 

.2231^-01 

24 

.25500+00 

,24^3-'+00 

74 

.22300-01 

-2164C-01 

25 

.22000+00 

.22063+00 

75 

.21600-01 

.21330-01 

26 

.20100+0° 

..I0793+OO 

76 

.21200-01 

,20730-01 

27 

.17700+00 

.17623+00 

77 

.20700-01 

.20233-01 

28 

.16000+00 

■  .15302+0° 

76 

.20300-01 

.la036-01 

29 

.14000+Q0 

.13.'.l2r'+0  0 

70 

.19900-01 

.19407-01 

30 

.12400+00 

.12246+00 

80 

.19500-01 

.lioQO-01 

31 

.11000+00 

.10826+00 

81 

.19200-01 

.18691-01 

32 

.10100+00 

.99984-01 

82 

.19000-01 

.18521-01 

33 

.95000-01 

.^3872-01 

83 

.18300-°! 

.18377-01 

34 

.890  0  0-01 

.88021-01 

P4 

.13500-01 

.1815c>-01 

35 

.32000-01 

.30444-01 

85 

.18300-01 

.18125-01 

3b 

•77500-°1 

.76642-01 

86 

.I8IOC-11 

.17839-01 

37 

.73500-01 

.72082-01 

87 

.17900-01 

,17872-01 

3° 

.70000-0! 

.69315-01 

88 

.17700-01 

.17,11-01 

39 

.66500-01 

.65177-01 

.   89 

.1750  0-01   ■ 

.17711-01 

4  0 

.63000-01 

.62344-01 

41 

.60500-01 

.5932.0-0.1 

90 
91 

.17400-01 

.17200-r 1 

17776-01 
,17007-01 

42 

.58000-01 

.57343-01 

9.2 

.17100-01 

17;i37-01 

43 

.5650  0-01 

.5550^-01 

93 

.17000-01 

1° 133-01 

4  4 

.'54500-0! 

.538  33-01 

94 

.16900-11 

,18216-01 

45 

•52000-ji 

.51002-01 

05 

.16800-01 

.18418-01 

46 

.49503-01 

.48799-01 

96 

.16700-01 

,  ip,o5^_03 

47 

.47000-01 

.46074-01 

97 

.16600-01 

.19310-01 

48 

.45000-01 

.44404-01 

98 

.16500-^1 

.19162-01 

49 

.43000-01 

.42071-01 

QFIH 

.16400-01 
F.MD  OF  RUN 

,13487-01 

MIFQ 


65 
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SPLCIFICATIONS 
IWTOTAL  RANGE  OF  DATAr=135 


NMf NUMBER  OF  DATA  REnri=128     IST*RT»STARTING  POTNT»=  i 


W   TS  CONSTANT 

ARBITRARY  UNITS 

INPUT  VECTOR   ZERQTH  ORDER   TRAPEZOID 


TRAPEZOTO+Pe'RT.  TP ^P^ZOTn+^PEPT.  NYSTROM  INTER, 


XKI) 


yyit; 


YKI) 


YK(I) 


Yi.m 


YM(T) 


7 

.607600+04 

.141759+03 

1800^2+0-3 

. 1 37279+03 

.116902+03 

.000000 

6 

.577700+0* 

.121642+03 

103465+03 

.  l??.2clc,+07 

.  1?5?03+03 

.131.1  23  +  03 

9 

.557700+04 

.112883+03 

139319+03 

.  l?47"?+03 

.119938+03 

.116242+03 

10 

.536100+04 

.967098+0? 

3594  75+0? 

.°A51^6+0? 

.107257+03 

,loao67+03 

11 

.522200+04 

.979292+02 

1 07472+03 

. 100704+0^ 

.989904+02 

,970185+02 

12 

.51.0900  +  04 

.915276+02      . 

883-O63+0? 

•°?6144+02 

.936497+02 

,9aa?0^+02 

13 

.493000+04 

.860158+02 

946683+02 

,o(i737(j+02 

.3°6208+02 

.8«6c;87  +  02 

11 

,452600+04 

.331039+02 

773628+0? 

,8?1603+02 

.833034+02 

,P44'»16+02 

15 

.469800+04 

.799909+02 

363451+0? 

. ^^9800+0? 

.824256+02 

.814134+02 

16 

.459400+04 

.782903+02 

711363+0? 

.764727+0? 

.78^733+02 

.7917?5+02 

17 

,4noooo+04 

,767<i37+02 

854438+0? 

.003  5o?+o? 

.784069+02 

.776704+02 

lb 

.433400+04 

.752085+02 

680«36+0? 

.733310+02 

.750852+02 

.759470+02 

19 

.  4?.8800+04 

.734904+0? 

323734+02 

.770241+02 

.751095+0? 

.7371 16+02 

20 

.418800+04 

.712085+0? 

646074+02 

.607680+02 

,71«654+0? 

.717^61+02 

21 

.414800+04 

.762621+02 

778096+0? 

.750934+02 

.745773+02 

.7^2880+02 

22 

.411(100  +  04 

,609377+02      . 

7)4  7146+0? 

.773215+02 

.783779+02 

.770O9o+02 

C  w 

.40530  0+04 

.894163  +  0?. 

872607+0? 

.PS8832+02 

.857960+0? 

.8H6668+02 

24 

.307900+04 

.100717+03 

915719+0? 

.°390  1.8+0? 

.°43?21+0? 

.9RRQQq+pj> 

25 

.300700+04 

» 111732+03 

109T62+0  3 

.107437+03 

.106197+03 

.107830+03 

2o 

.346800+04 

.115263+03 

113601+03 

.  1.1  3532+03 

.ll?69?+03 

.116810+03 

27 

.294700+04 

.111571+03 

1.16^26+03 

.1  I.4507+03 

.ll?17?+03 

.1168?o+03 

2c 

.222100+04 

.056671+02 

1  062.1  5+03 

.101151+03 

.977818+02 

,985?68+0? 

29 

.13670  0+04 

.445148+02 

651186+02 

.651002+0? 

.6491 15+02 

.622325+02 

30 

.791000+03 

.153256+0? 

239110+02 

.270171+02 

.309654+02 

.275775+02 

31 

•5ft3r00+03 

.500530+01 

674 r 2  8+ 01 

.°02374+01 

.113056+02 

.921 071+01 

3£ 

.558000+03 

,318811+01 

327031+01 

.3821 '4+ 01 

.4673Ro+01 

.380135+01 

33 

.577000+03 

,418990+01 

310590+0] 

.349464+03 

.3°.1836+01 

.363?08+01 

34 

.609000+03 

.586739+01 

527389+01 

.511039+01 

.500126+01 

.506047+01 

35 

.628000+03 

.750582+01 

646089+01 

.661137+01 

.66^931+01 

•677?17+01 

36 

.639000+03 

.896337+01 

855P75+0.1 

,O339°«+01 

.821320+01 

.8?7R3U+01 

37 

.630000+03 

.958440+01 

937600+ni 

.c'307n?  +  n1 

.o?6041+01 

,9?«786+ol 

3o 

.614000+03 

.104151 +0? 

979279+01 

.003075+QI 

.006^81+01 

.100^53+02 

39 

.589000+03 

.108794+02 

110374+n? 

.1  n7b'<-0  +  02 

.106464+02 

.107^48+02 

40 

.551000+03 

.111942+02 

1^64 16+0? 

.  1 08917+02 

.100171+02 

.111611+0? 

41 

.500000+03 

.105743+02 

1 17468+C? 

.11 1718+02 

.108808+0? 

. 10°6R4+02 

42 

.43300  0+03 

.930339+01 

940185+01 

.075986  +  0.1 

.op?5«4.  +  oi 

,9O?-«00  +  01 

43 

.368000+03 

.772316+01 

9?0492+01 

.074302+01 

.856978+01 

,8r'?07l  +  01 

4  4 

.305000+03 

.608601+  )1 

624.140+01 

.£.£,83c?+0! 

.683034+01 

.630561+01 

45 

.252000+03 

.449584+01 

5Q3062+0  5 

.560415+01 

.536894+01 

,5?474l+ni 

46 

.207000+03 

.298367+01 

3061 06+01 

.361353  +  0.'! 

.36Q00?+01 

.366807+01 

67 


47 

.  l"7Mino  +  Q3 

.185556+01 

.290628+01 

.2581*4+01 

.251962+01 

.235690+01 

48 

.'15&C00+03 

.11701 4+01 

.R04843+C0 

.127685+01 

.147460+01 

.145925+01 

49 

.150000+03 

.769560+00 

•  153543+n 

.115838+01 

.106458+01 

.924647+00 

5C 

.l4<~O00  +  03 

.658164+00 

.368735-0? 

.477137+00 

.665146+00 

.6°962l+00 

51 

.152000+03 

.801045+00 

.131264+01 

.923950+00 

.806459+00 

.732213+00 

52 

.165000+03 

.906958+00 

.269449+00 

.665817+00 

.738759+00 

.842497+00 

53 

.158000+03 

.100458+01 

.152447+01 

.1.14534+01 

.102107+01 

.937598+00 

54 

.160000+03 

.114847+01 

.484689+00 

.879247+00 

.102744+01 

.109*59+01 

55 

.165000+03 

.154620+01 

.181225+01 

.150231+01 

.138592+01 

.136143+01 

56 

.163000+03 

.145475+01 

.128015+01 

.142703+01 

.146365+01 

.149111+01 

57 

.160000+03 

.163028+01 

.162935+01 

.1571^6+01 

.155924+01 

.154507+01 

58 

.159000+03 

.166619+01 

.163120+01 

.164255+01 

.164419+01 

.164009+01 

59 

.154000+03 

.180312+11 

.170117+01 

.172349+01 

.17^044+01 

.173020+01 

60 

.150000+03 

.183011+01 

•1°0606+01 ' 

.1*4610+01 

.182246+01 

.1*2221+01 

61 

.142000+03 

.189638+01 

.175516+01 

.182688+01 

.184750+01 

.187047+01 

62 

.134000+03 

.1863C6+01 

.203561+01 

.193135+01 

.188777+01 

.190305+01 

63 

.122000+03 

.176046+01 

.169050+01 

.177434+01 

.178822+01 

.1*3029+01 

64 

.109000+03 

.148404+01 

.184842+01 

.170064+01 

.164033+01 

.162052+01 

65 

.949000+02 

.119156+01 

.111966+01 

.126509+01 

.131772+01 

.133342+01 

66 

.849000+02 

.981093+00 

.126346+01 

.114537+01 

.110994+01 

.108476+01 

67 

.758000+02 

.759434+00 

.698729+00 

.813085+00 

.851320+00 

.862485+00 

66 

.682000+02 

.551062+00 

.820139+00 

.71.0212+00 

.679618+00 

.651927+00 

69 

.635000+02 

.438280  »-00 

.231986+00 

.423776+00 

.473783+00 

.4*658*+00 

70 

.598000+02 

.279286+90 

.594575+on 

.437380+00 

.388945+00 

.345623+00 

71 

.581000+02 

.237834+00 

-.360034-01 

.160372+00 

.234530+00 

.253"*76+00 

72 

.537000+02 

.235595+00 

.5il672+00 

.328367+00 

.271787+00 

.237143+00 

73 

.595000+02 

.275524+00 

-.404817-01 

.ittR-'V+OO 

.221258+00 

.250751+00 

74 

.599000+02 

,247940+00 

.591630+00 

.371665+00 

.302165+00 

.261446+00 

75 

.612000+02 

.356068+00 

-.956501-01 

.169453+00 

.2c7e03+00 

.300740+00 

76 

.622000+02 

.328169+00 

.*P 7785+0  0 

.497341+00 

.393691+00 

.340124+00 

77 

.620000+02 

.433273+00 

-.151446+00 

.203332+00 

.324941+00 

.3*7075+00 

76 

•637COC+02 

.465672+00 

.101799+01 

.63897Q+00 

.505755+00 

,445721+00 

79 

.609000+02 

.446827+00 

-.866470-01 

.2752*4+00 

.390367+00 

.4616R2+00 

80 

.617000+02 

.541131+00 

.980301+00 

.656086+00 

.545470+00 

.400597+00 

81 

.590000+02 

.476466+00 

.  101of,2+'"  0 

."73103+00 

.459122+00 

.505120+00 

82 

.569000+02 

.490132+00 

.851Q10+00 

.605876+00 

.525614+00 

.479348+00 

83 

.550000+02 

.451570+00 

. 129254+00 

.3569*5+00 

.436188+00 

,4"'4625  +  00 

84 

.535000+02 

.524485+00 

.773885+00 

.583313+00 

.51*267+00 

.493457+00 

85 

.512000+02 

.458518+00 

•275084+^0 

.419363+00 

.^62768+00 

.489053+00 

86 

.475000+02 

.  4147053+rin 

.641952+00 

.515841+00 

.475609+00 

.464204+00 

87 

.4^5000+02 

.413575+00 

.252153+00 

.370927+00 

.410215+00 

.433371+00 

86 

.423000+02 

.396640+00 

.574997+00 

.461738+00 

.421547+00 

.404407+00 

89 

.394000+02 

.319288+00 

.218283+00 

.311404+00 

.342740+00 

.357776+00 

90 

.366000+02 

.307685+00 

.420293+00 

.3490*9+00 

.326109+00 

.315393+00 

91 

.346000+02 

.243915+00 

.195077+00 

.248892+00 

.265388+00 

.273774+00 

92 

.320000+02 

.206338+90 

.^°2763+00 

.247669+00 

.233550+00 

.221680+00 

93 

.308000+02 

.158938+00 

.119924+0  0 

.161733+00 

.178021+00 

.181672+00 

94 

.2^8000+02 

.163684+00 

. l«7O52+r0 

.1735^+00 

.166597+00 

.l61°59+00 

95 

.295000+02 

.138160+00 

.129416+0" 

.143754+00 

.147630+00 

.14*271+00 

9b 

.286000+02 

.126659+00 

,146°04+00 

.1.37241  +  00 

.135c;6c;+00 

.13113^+00 

68 


97 

.2^:6000  +  02 

,12Pc»47+00 

.106414+00 

.120674+00 

.126791+00 

.l?0533+00 

9H 

2£ 60 00+02 

.  1.42012+00 

.151480+00 

.140813+00 

.136082+00 

.135279+00 

99 

2rtS000  +  02 

,123136+00 

.13254 4+ CO 

.132564+00 

.132417+00 

.130"'77+00 

1 0  0 

2.35000  +  02 

.  13344 0  +  0n 

.11 3727+O0 

.133434+00 

.127932+00 

.127739+00 

101 

289000+02 

,137271+30 

.153154+00 

.141239+00 

.137673+00 

.135051+00 

102 

>2,J0000  +  02 

.153683+00 

.121389+00 

.137443+00 

.143121+00 

.145601+00 

103 

2'-)300n  +  02 

,lb32(J2  +  0n 

.185977+00 

.167647+00 

.161362+00 

.157793+00 

104 

243000+02 

.176538+00 

.  I.aa583+'t0 

.160136+00 

.166733+00 

,16Q258+00 

105 

.293000+02 

.188575+Pfj 

.212488+00 

.192534+00 

.186941+00 

.134700+00 

|06 

2^4000+02 

.220392+00 

.164661+00 

.19137^+00 

.193591+00 

.203830+00 

107 

287000+02 

.202513+00 

.277122+00 

.23351  1.  +  00 

.21°097+00 

.210248+00 

10a 

2^1000+02 

.237065+00 

.127914+00 

.189186+00 

.210914+00 

.22.1963+00 

109 

278000+02 

.242651+00 

.346215+00 

.275310+0" 

.250024+00 

.241802+00 

110 

261000+02 

.2474  81+0  0 

. 130086+PP 

.209739+00 

.232567+00 

,248°23+00 

111 

249000+02 

,240043+00 

.355375+00 

.281135+00 

.253574+00 

.244534+00 

112 

2  24000+02 

,197238+oq 

.1.24221  +  00 

. 1 87160+00 

.207392+00 

.217560+00 

133 

2 070  0 0+02 

.176144+00 

.2702^5+00 

.214546+00 

.196812+00 

.186424+00 

114 

1  ^2000  +  02 

.148384+00 

.320339-01 

•1355°l+00 

.153449+00 

.161655+00 

13b 

17900H+02 

,12°0B'l  +  O0 

.214735+00 

.164063+00 

.147731+00 

.133613+00 

116 

170000+02 

.111267+00 

.434341-01 

.On5953-0l 

.111635+00 

.11°?00+00 

117 

162000+G2 

.926072-01 

.179101+10 

.127658+00 

.111084+00 

.100585+00 

118 

1^6000+02 

.Bb0^84-01 

.61175&-02 

.612730-01 

.808624-01 

.878161-01 

119 

1P8000+02 

.844897-01 

.164083+00 

.111224+00 

.940922-01 

,83000°-01 

120 

lS«000+02 

.857636-f1] 

.489619-0* 

.5A3849-01 

.772763-01 

.842332-01 

121 

.162000+02 

.  104273+00 

.166641+0  0 

.1 18894+00 

.103765+00 

. 961^77-01 

122 

.164000+02 

.119697+00 

.4  19046-01 

.086250-01 

.103699+00 

.112752+00 

123 

164  000+02 

.129221+00 

. iQ7aP9+oo 

.148302+00 

.131465+00 

.122800+00 

124 

.159000+02 

.12.5052+00 

.609521*01 

.105075+00 

.121152+00 

.128343+00 

12b 

lbOCOO+02 

,159223+00 

.169151+00 

.157809+00 

.146820+00 

,14370Q+00 

126 

154000+02 

.145321+00 

.1292«5+00 

.144613+00 

.148344+00 

.151653+00 

127 

lui+QOO  +  02 

.154748+np 

.161346+00 

.153605+00 

.181714+00 

.150595+00 

128 

.139000+02 

.148464+00 

. 148150+00 

,15045/i  +  QO 

,151080+00 

.153043+00 

129 

1260CO+02 

.  Ib5336  t-CO 

•  148770+00 

.150860+00 

.151.849+00 

.152258+00 

130 

.122000+02 

.154366+00 

.1618,03+0  0 

.157198+00 

.155001+00 

.14283°+00 

131 

109000+02 

.149861+00 

.1.46839+00 

.150355+00 

,152900+00 

.131625+00 

132 

103000+02 

.173593+pn 

.152083+00 

.158779+00 

.161875+00 

.191206+00 

133 

976000+01 

.18Q43r,+nfj 

.194302+00 

.185776+00 

.150314+00 

.370020+00 

134 

951000+01 

.184567+00 

.1^4567+^0 

.neoooo 

.ooooop 

.562712+00 

69 


BACK  CHECKS 
TMPUT      RECOMPUTED 


XT(I) 


ZKT) 


7  - 

.60760+04 

.60833+04 

57 

.16000+03 

.16285+03 

ft 

.57770+04 

.57555+04 

58 

.  1 5900  +  0  "^ 

.15931+03 

9 

.55770+04 

.55991+04 

59 

.lS'+oo+o's 

.15^24+03 

10 

.53610+04 

.53151+04 

60 

.15000+13 

,15n??+03 

11 

.52220+04 

.52649+04 

61 

.14200+03 

.1425^+03 

12 

.51090+04 

.50973+04 

62 

.13400+03 

. 1340°+O3 

13 

,<t9300  +  04 

.49539+04 

63 

.12200+m 

.12225+03 

14 

.48260+04 

,4ft071+04 

64 

.10900+03 

.10986+03 

15 

.469S0+04 

./j700^  +  04 

65 

.9490  0+0? 

."5175+02 

16 

.45940  +  0'* 

. 4573^+04 

66 

.84900+02  - 

.r^os^  +  O? 

17 

.44900+04 

.44081 +04 

67 

.75300+0? 

.7541^+02 

la 

.43840+04 

.4373?  +  0<J 

68 

.68200+0° 

.690-(5  +  C2 

19 

.42980+ J4 

.42^83+04 

69 

.63500+0° 

.63201+02 

?o 

.41680+04 

.41798+04 

70 

.59800+02 

.60210+02 

?i 

.41480+04 

.41841+04 

71 

.58100+0° 

.58481+02 

22 

.41140+04 

•410&4+04 

75 
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Figure  1.  Scattering  from  melt-crystallized  polyethylene,  at  very  low 
angles  (20^0.008  radians).  Weighting  function  is  constant. 
Solid  line-Nystrom  (forward)  interpolation. 
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Bottom  curve    (represented   by   Q     points)    give   the   input 
data. 

O    -   Nystrom    (forward)    interpolation. 

X   -   rectangle  rule  approximation. 
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Figure  2.   Data  from  the  same  experiment  as  in  Figure  1,  but  for 

scattering  angles  in  the  range  of  (0.008<26<:.02  radians) 
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O  -  rectangle  rule  approximation. 


X  -  trapezoid  rule  approximation,  corrected  by  smoothing  the 
results. 
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Figure  3,   Unsmearing  of  experimental  data  which  were  not  smoothed  out 
for  random  noise  effects.   Bottom  curve  (represented  by  □ 
points)  show  the  unsmoothed  input  data.   Solid  line-Nystrom 
(forward)  interpolation. 
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O  -  rectangle  rule  approximation. 

X  -  trapezoid  rule  approximation,  corrected  by  smoothing 
the  results. 
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Figure  4.   Scattering  from  melt-crystallized  polyethylene.   N=1UU, 
trapezoid  weighting  function.   Bottom  curve  (represented 
by    points)  show  the  input  data.   Solid  line-Nystrom 
(forward)  interpolation. 
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Figure  5.   Scattering  from  triblock  copolymer.   Constant  weighting 
function.   Bottom  curve  shows  the  input  data.   Upper 
curve  shows  desmeared  data,  Nystrom  (forward)  inter- 
polation. 
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Figure  6.  Scattering  intensities  by  spheres  with  uniform  electron 
densities.  Solid  line-exact  results.  Points  represent 
computed  results,  based  on  Nystrom  (forward)  interpolation, 
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